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ABSTRACT 


\ 

--■'A  computational  method  coupling  the  three  degrees  of  freedom  flight 
mechanics  equations  and  the  two  dimensional  Navier-Stokes  equations  was 
developed  which  could  be  used  to  predict  the  flight  path  of  a  free 
falling,  autorotating,  two  dimensional  flat  plate.  The  two  dimensional 
incompressible  Navier-Stokes  equations  were  cast  in  a  body  fixed 
coordinate  system.  The  corresponding  velocities  were  cast  in  an 
inertial  reference  system.  The  equations  were  represented  by 
backward-time-central -space  finite  differences  and  solved  using  a 
successive-over-relaxation  iteration  technique.  The  resulting 
aerodynamic  coefficients  were  entered  into  the  three  degrees  of  freedom 
flight  mechanics  equations.  The  system  of  ordinary  differential 
equations  was  solved  using  an  Adams  open  formula  to  predict  the  movement 
of  the  plate.  New  boundary  conditions  for  the  Navier-Stokes  equations 
solver  were  derived  from  the  movement  of  the  plate.  ■  The  process  was 
repeated  to  advance  the  solution  in  time. 

The  computational  method  was  used  to  calculate  the  flew  field 
around  a  flat  plate  forced  to  rotate  at  nondimensional  angular 
velocities  of  1.0,  2.0,  and  4.0.  ^A  calculation  was  performed  for  a 
plate  with  one  degree  of  freedom  (^rotation)  which  simulated  a  plate 
mounted  on  a  frictionless  bearing  in  a  wind  tunnel.  The  plate  was  found 

to  undergo  pinned  autohotation  at  a  mean  nondimensional  rotational 

/ 

velocity  of  1.58.  —^Finally,  the  flight  path  of  a  free  falling 
autorotating  plate  was  predicted  using  the  computational  procedure.  The 
validity  cf  the  overall  approach  was  demonstrated  by  comparison  with 
experiment. 


SECTION  I 


INTRODUCTION 


The  USAF  has  lost  hunareds  of  millions  of  dollars  worth  of  aircraft 
to  accidents  which  were  attributed  to  aircraft  departure  during  stall 
which  resulted  in  a  spin  (Ref  1).  Examination  of  the  accident 
statistics  for  a  large  number  of  NATO  combat  aircraft  showed  that 
approximately  20  percent  of  the  aircraft  losses  in  combat  training  were 
caused  by  stall  departure  and  spinning  accidents  (Ref  2).  Many  of  the 
accidents  may  have  been  prevented  if  the  spin  characteristics  of  the 
aircraft  had  been  better  understood.  However,  predicting  the 
characteristics  of  a  spinning  aircraft  is  very  involved.  Conventional 
wind  tunnel  testing  of  aircraft  models  does  not  accurately  simulate 
spin.  Experience  has  indicated  aircraft  spin  can  be  investigated  safely 
and  at  a  comparatively  moderate  cost  by  testing  small  dynamic  models  in 
a  vertical  wind  tunnel  (Ref  3).  However,  because  of  the  many  variables 
ir>  a  spin,  interpretation  of  the  results  for  application  to  a  corre¬ 
sponding  airplane  is  difficult  and  may  only  show  general  trends. 
Ultimately,  it  is  necessary  to  turn  to  full  scale  flight  testing  which 
is  extremely  expensive.  One  is  inclined  to  turn  to  numerical  techniques 
to  analyze  aircraft  spin  because  the  per  unit  computing  costs  are 
decreasing  rather  than  increasing.  Unfortunately,  the  numerical 
aerodynamic  analysis  of  an  aircraft  during  spin  is  extremely  complex. 

The  flight  condition  of  a  "spinning"  aircraft  was  defined  by 
Bairstow  (Ref  4)  as  a  steady  spiral  descent  with  the  mean  angle  of 
attack  of  the  wing  greater  than  stalling.  This  latter  element  infers 
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ilterf  is  unsteady,  separated  flow  about  the  aircraft  which  makes  it 


difficult  tc  mathematically  predict  the  spir  characteristics  of  the 
aircraft.  Quasi -steady  techniques  developed  tc  predict  aircraft  spir, 
are  incapable  of  correctly  defining  the  true  dynamics  of  the  flow  field 
(Ref  5).  Only  solutions  of  the  Navier-Stokes  equations  can  accurately 
predict  unsteady,  separated  flow  fields. 

To  mathematically  investigate  aircraft  spin,  a  six  degrees  of 
freedom  trajectory  program  should  be  combined  with  a  three  dimensional 
Navier-Stokes  program.  The  large  computer  memory  requirements  and  long 
run  times  make  this  task  impractical  on  today's  computers  (Ref  6). 
While  one  awaits  larger  computers  to  precisely  calculate  aircraft  spin, 
one  can  search  for  a  simpler  problem  which  can  be  accurately  solved  on 
present  computers.  Developing  and  demonstrating  a  procedure  to  solve  a 
simpler  problem  which  incorporates  the  essential  elements  of  a  spinning 
aircraft  will  be  a  step  towards  solving  the  complete  aircraft  spin 
problem. 

The  free  falling,  autorotating,  flat  plate  incorporates  the 
essential  elements  (unsteady  aerodynamics,  flight  mechanics,  and 
autorotation)  of  a  spinning  aircraft.  When  a  rectangular  piece  of  paper 
(e.g.  a  computer  card)  is  dropped,  it  falls  obliquely  as  shown  in  Figure 
la  and  not  vertically  as  might  be  expected.  The  flow  field  around  the 
free  falling  flat  plate  is  unsteady  and  separated.  Based  on 
experimental  observations,  the  autorotating  flat  plate  can  be  modeled  in 
two  dimensions  (Ref  7).  The  author  feels  the  flight  path  of  a  free 
falling,  autorotating,  flat  plate  can  be  numerically  predicted  using 
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present  computers. 


The  purpose  of  this  research  was  to  develop  and  apply  a  method  for 
solving  the  three  degrees  of  freedom  flight  mechanics  equations.  These 
equations  could  be  used  to  predict  the  flight  path  of  a  free  falling, 
autorotating,  two  dimensional  object.  The  three  degrees  of  freedom  are 
vertical  translation,  horizontal  translation,  and  rotation.  The 
unsteady,  nonlinear  aerodynamic  forcing  functions  (lift,  drag,  and 
moment)  in  the  flight  mechanics  equations  were  simulaneously  calculated 
by  solving  the  two  dimensional  incompressible  Navier-Stokes  equations. 
Past  attempts  failed  to  predict  the  flight  path  of  a  free  falling, 
autorotating,  flat  plate  because  the  aerodynamic  forces  were  linearized 
and  assumed  to  be  quasi-steady  (Ref  7).  Jhe  method  of  solution  was 
validated  by  predicting  the  flight  path  of  a  free  falling  flat  plate  and 
comparing  it  to  experimental  data.  To  the  author's  knowledge,  this  was 
the  first  attempt  to  use  an  unsteady,  nonlinear  aerodynamic  prediction 
method  to  calculate  the  aerodynamic  forces  included  in  the  three  degrees 
of  freedom  flight  mechanics  equations.  This  research  will  assist  in  the 
future  analysis  of  the  complete  aircraft  spin  problem. 
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SECTION  II 


BACKGROUND 


This  section  briefly  outlines  the  history  of  the  phenomenon  of  the 
free  falling,  autorotating,  flat  plate.  The  phenomenon  was  first 
described  in  the  middle  of  the  last  century,  but  early  investigators 
were  unable  to  explain  the  autorotating  plate  for  two  primary  reasons. 
First,  detailed  experimental  data,  being  difficult  to  obtain,  were  not 
available.  Secondly,  scientists  attempted  to  calculate  the  flow  field 
around  the  flat  plate  using  linear,  quasi-steady  techniques.  This  led 
to  erroneous  conclusions  because  the  autorotating  phenomenon  was  caused 
by  the  neglected  unsteady,  nonlinear  effects.  Fortunately,  in  the  past 
few  years,  detailed  experimental  data  and  accurate  nonlinear,  unsteady 
aerodynamic  prediction  techniques  have  aided  in  satisfactorily 
explaining  autorotation  of  free  falling  flat  plates. 

Maxwell  (Ref  8),  in  a  brief  theoretical  study  written-  in  1854,  was 
probably  the  first  person  to  attempt  to  describe  the  phenomenon  of  an 
autorotating  rectangular  plate  falling  through  the  air.  Maxwell's 
explanation  of  autorotation,  although  inaccurate,  is  summarized  with  the 
aid  of  Figure  la.  At  position  1,  the  plate  is  parallel  to  the  flight 
path  which  corresponds  to  the  least  drag.  The  plate  has  its  maximum 
drag  at  postion  3  because  it  is  normal  to  the  flight  path.  The  plate 
therefore  has  a  higher  velocity  at  position  2  than  at  position  4.  Since 
a  plate  always  tends  to  place  itself  normal  to  the  flight  path,  the 
torque  at  position  2  Is  in  the  direction  of  rotation  and  the  torque  at 
postion  4  is  opposite  to  the  direction  of  rotation.  Hence,  the  torque 
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at  postion  2  v.t.ich  supports  rotation  is  larger  than  the  adverse  torque 
at  postior,  4  which  resists  rotation.  Even  though  Maxwell's  explanation 
of  the  phenomenon  was  inaccurate,  several  French  and  German  papers 
written  between  1880  and  1900  essentially  restated  it  (Ref  7).  They 
assumed  the  total  torque  could  be  divided  into  a  quasi-steady  part  which 
always  acted  in  the  direction  of  rotation  and  a  quasi-steady  part  which 
always  acted  opposite  to  the  direction  of  rotation. 

It  was  not  until  around  1905  that  Rianbouchinskiy  realized 
Maxwell's  quasi-steady  theory  would  not  account  for  a  flat  plate 
autorotating  about  a  fixed  axis  (Ref  7).  This  may  have  suggested  to  him 
Maxwell's  theory  was  not  valid  for  the  free  falling  plate  either.  He 
recognized  Maxwell's  assumption  that,  in  position  1,  the  movement  of  the 
flat  plate  was  faster  than  it  was  in  position  3,  which  was  not  valid  for 
a  flat  plate  autorotating  about  a  fixed  axis..  Riabouchinskiy  iRef.  9) 
explained,  that  in  position  1,  the  streamlines  at  the  retreating  edge  of 
the  plate  are  more  curved  than  in  position  3,  thus,  due  to  the  higher 
suction  effect,  the  torque  favorable  to  rotation  was  greater.  It  was 
not  until  several  years  later  that  he  realized  the  moment  of  inertia 
must  be  large  enough  to  overcome  the  period  of  adverse  torque  (Ref  7). 

In  1941,  Dupleich  (Ref  10)  published  the  first  experimentally 
obtained  trajectories  of  free  falling  plates  in  sir  and  water. 
Dupleich's  explanation  of  autorotation  was  adopted  from  Maxwell  and 
depended  on  the  incorrect  idea  that  the  influence  of  rotation  on  the 
total  torque  was  always  opposite  to  the  direction  of  rotation. 


F.K.C.  Smith  (Ref  11)  ir  1953  put’ished  the  first  accurate 
explanation  of  the  free  foiling  flat  plate,  he  realized  the  asymmetry 
cue  to  rotation  was  the  key  tc  explairinc  the  autcrotation  and  this 
asymmetry  was  caused  by  hysteresis  during  acceleration.  Smith  explained 
that  at  position  2,  as  shown  in  Figure  lb,  the  plate  was  at  a  large 
angle  of  attack,  producing  a  force  R.  This  force  R  causes  a  nose  up 
moment  which  tended  to  increase  the  angle  of  attack.  Rotation  was 
rapid;  consequently  the  boundary  layer  was  thinner  than  it  would  have 
been  for  the  steady  state  condition  and  the  stall  was  delayed  to  high 
values  of  lift  coefficient.  At  position  3,  the  plate  was  probably 
completely  stalled.  At  position  4,  it  continued  to  be  stalled. 
However,  between  positions  4  and  5,  attached  flow  over  the  trailing  edge 
was  established.  Reattachment  was  delayed  _to  angles  below  the  static 
stall  angle  because  the  thick  wake  was  a  retarding  influence.  Hence  the 
downward  lift  impulse  was  less  than  the  ypward  lift  impulse.  -The 
clockwise  rotation  impulse,  and  the  inertia  of  the  plate,  carried  it 
past  the  region  of  retarding  impulse.  The  motion  from-  5  to  9  was 
essentially  a  repetition  of  that  from  1  to  5. 

In  the  Tate  1960 ’ s  Bustamante  and  Stone  (Ref  12)  experimentally 
investigated  autorotation  of  flat  plates  at  subsonic  and  hypersonic 
speeds.  They  found,  although  not  conclusively,  that  plates  could 
autorotate  in  hypersonic  flow.  They  suggested  that  autorotation  of 
wings  with  a  symmetric  cross  section  is  due  to  a  large  vortex  shed  from 
the  retreating  face  of  the  wing.  Iversen  (Ref  13)  was  able  to  show  the 
same  result  theoretically  in  1969. 


In  IS7C,  E.H.  Smith  (Ref  14)  published  the  first  detailed 
experimental  data  on  the  lift.,  drag,  moment,  angular  acceleration,  and 
angular  velocity  as  a  function  of  time  for  a  flat  plate  autorotat.ing 
about  a  fixed  axis  in  a  wind  tunnel.  Smith  studied  the  corresponding 
flow  patterns  using  smoke  and  tufts.  He  also  experimentally 
investigated  free  falling  wings.  He  found,  for  Reynolds  numbers  above 
4000,  the  average  lift  and  drag  coefficients  were  comparable  to  those 
observed  in  the  fixed  axis  tests. 

In  1974,  Kry  and  List  (Ref  15)  experimentally  determined  the  range 
of  autorotation  in  which  a  quasi-steady  technique  could  be  applied  to 
rotating  oblate  spheroids.  They  found,  if  the  Reynolds  number  was  less 
than  150,000,  the  moment  predicted  for  the  oblate  spheroid  was  valid  for 
rotation  rates  about  the  major  axis  equivalent  to  Strouhal  numbers  less 
than  0.04.  Under  their  conditions,  the  moment  acting  on  an  obljte 
spheroid  in  a  steady  flow  differed  by  less  than  10%  from  the  moment 
acting  on  a  free  spheroid  with  the  same  instantaneous  orientation.  List 
actually  computed  rates  of  autorotation  by  means  of  the  quasi -steady 
approach. 

In  1978,  Lugt  (Ref  7)  numerically  predicted  the  autorotation  of  a 
thin  elliptic  cylinder  rotating  about  a  fixed  axis.  Lugt  transformed 
the  elliptic  cylinder  from  the  physical  plane  to  a  computational  plane 
where  he  solved  the  Navier-Stokes  equations.  The  equations  which  were 
valid  for  laminar,  incompressible,  two  dimensional  flow  were  formulated 
in  terms  of  vorticity  and  the  stream  function.  The  initial  conditions 
consisted  of  the  potential  flow  solution  and  a  vorticity  sheet 
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distributee  alone  the  surface  of  the  cylinder  to  rnoael  a  no  slip 
condition.  Loot  significantly  reduced  the  computer  time  necessary  to 
arrive  at  c  solution  by  assuming  the  cylinder  autoroteted  at  a  constant 
angular  velocity.  He  predicted  the  autorotation  of  the  elliptic 
cylinder  and  showed  autorotation  was  caused  by  a  complicated  interplay 
of  vortex  shedding,  boundary  layer  hysteresis,  and  vorticity  generation 
around  the  edges.  Autorotation  and  the  aerodynamic  forces  on  the 
ellipse  predicted  by  Lugt  correlated  very  well  with  the  high  Reynolds 
number  and  three  dimensional  experimental  data  obtained  by  E.H.  Smith 
(Ref  14). 

Researchers  attempted  to  explain  the  free  falling,  autorotating. 


flat  plate  as  early  as  1854.  It  was  not  until  1953  that  A.M.O.  Smith 

4 

was  able  to  explain  the  phenomenon.  To  the  author's  knowledge,  no  one 
has  been  able  to  accurately  predict  the  flight  path  of  a  free  falling, 
autorotating,  flat  plate.  In  this  study  a  computational  method  was 
developed  which  coupled  the  three  degrees  of  freedom  flight  mechanics 

equations  with  the  two  dimensional  Navier-Stokes  equations.  This  j 

i 

numerical  method  was  used  to  predict  the  flight  path  of  a  free  falling, 

I 

autorotating,  flat  plate.  | 
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SECTION  III 


NAVIER-STOKES  ECUATIOKS  SOLVER 


A  numerical  two  dimensional  incompressible  Navier-Stokes  equations 
solver  was  developed  to  calculate  the  aerodynamic  characteristics  of  a 
flat  plate.  The  flow  field  was  modeled  by  a  finite  number  of  points. 
The  finite  difference  representation  of  the  Navier-Stokes  eouetions  was 
solved  at  each  of  these  points  to  obtain  an  x  velocity  component,  a  y 
velocity  component,  and  a  static  pressure.  The  aerodynamic 
characteristics  of  the  plate  were  obtained  by  integrating  the  forces 
around  the  plate. 

A.  Grid  Description 

The  aerodynamic  forces  acting  on  the  plate  were  obtained  by  solving 
the  Navier-Stokes  equations  using  a  finite  difference  technique..  This 
required  the  flow  field  to  be  modeled  by  discrete  points.  An 
incompressible  flow  field  should  have  an  infinite  outer  boundary  because 
each  point  in  the  flow  field  influences  every  other  point  in  the  flow 
field.  An  infinite  outer  boundary  is  not  feasible  nor  is  it  necessary 
in  numerical  modeling  of  a  flow  field.  Previous  researchers  (Refs  16 
and  17)  found  placing  the  outer  boundary  10  chords  away  from  the  body 
closely  approximated  an  infinite  boundary.  Therefore  the  outer  boundary 
was  placed  at  10  chords  throughout  this  study. 

The  Navier-Stokes  equations  were  solved  over  a  rectangular  physical 
domain  in  this  study.  The  domain  was  modeled  by  I  MAX  points  distributed 
horizontally  in  the  x  direction  and  JMAX  points  distributed  vertically 
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it,  the  y  direction.  The  grid  point:  :r  the  physic: '  plane  were 
clustered  near  the  plate  v.'here  large  velocity  arc  pressure  gradients 
were  anticipated.  The  grid  point  clustering  was  arranged  so  it  did  r.ct 
favor  any  particular  flow  direction  because  the  freestrean  velocity 
could  originate  from  any  direction.  Ir,  addition,  the  points  were 
distributed  so  they  could  be  transformed  into  a  computational  domain 
with  uniform  grid  spacing  in  both  directions. 

Many  investigators  used  wrap  around  grid  systems  when  solving  the 
Navier-Stokes  equations  over  airfoils.  A  wrap  around  grid  system 
(0-Grid)  was  not  used  in  this  study  because  it  could  present  numerical 
difficulties  at  sharp  edges.  Hodge  (Ref  16)  found  the  large  leading 
edge  radius  on  an  airfoil  did  not  present  any  difficulty  whereas  the 
sharp  trailing  edge  was  a  constant  source  of  numerical  difficulty. 
Similar  difficulties  were  anticipated  on  the_ edges  of  the  flat  plate 
investigated  in  this  study. 

Hodge  (Ref  18)  found  an  exponentially  stretched  grid  inferior  to  an 
optimized  grid  in  solving  the  flow  field  around  an  airfoil.  The 
optimized  grid  minimized  the  truncation  error  assuming  steady  Blasius 
profiles.  In  addition,  Hodce  also  found  the  maximum  grid  spacing  near 
the  outer  boundary  necessary  to  obtain  acceptable  accuracy  was 
approximately  one  chord.  However,  combining  an  exponentially  stretched 
grid  with  a  very  small  spacing  near  the  plate,  and  a  maximum  grid 
spacing  or  approximately  one  chord  on  the  outer  boundary  required  many 
grid  points. 
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A  stretched  Cartesian  coordinate  syster.  was  selected  as  a 
compromise.  It  avoided  the  numerical  difficulties  of  a  wrap  crounc  gric 
system.  It  was  unnecessary  to  accurately  resolve  the  bounGary  layer  or 
the  plate  during  this  study  because  the  pressure  forces  produced  by  the 
large  eddy  structure  were  expected  to  dominate  the  problem  (Ref  19). 
Therefore  it  was  decided  to  distribute  the  grid  points  to  more 
accurately  resolve  the  large  eddy  structure.  An  adaptive  grid  is 
another  approach  to  optimizing  the  grid,  however,  that  is  beyond  the 
scope  of  this  study. 

The  Cartesian  grid  simplified  the  transformation  metrics  because 
each  variable  in  the  computational  domain  was  a  function  of  only  one 
variable  in  the  physical  plane  as  seen  in  equations  (1)  and  (2). 

X  =  X  (  o 

y  =  y(n) 

The  grid  was  generated  with  analytic  functions.  This  presented  the 
option  of  calculating  the  transformation  metrics  with  numerical  or 
analytical  derivatives.  Both  types  of  derivatives  were  investigated  and 
are  discussed  in  Section  VI. 

The  thickness  of  the  plate  was  used  to  determine  the  minimum  grid 
spacing  in  the  y  direction.  It  would  have  been  computationally 
prohibitive  to  accurately  model  this  mimimum  thickness.  The  minimum 
grid  spacing  was  selected  to  produce  a  grid  with  an  outer  spacing  of 
approximately  one  chord  and  a  total  number  of  points  which  permitted  the 


(1) 

(2) 
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profiler;  tc  be-  solvec  cr  e  CYBER  750  computer.  The  ectue’-  thickness  o' 
the  plate  was  ret  felt  tc  te  a  key  feature  cf  autcrctetior.. 


A  patched  exponential  stretching  scheme  of  the  following  form  was 
used  tc  smoothly  distribute  the  grid  points  in  the  physical  plane 


(e  y  ' -1) 
C 

(e  y-l) 


for  j  =  1,  JHi}-/: 


where  is  the  stretching  coefficient  in  the  y  direction  (Ref  20).  The 
remaining  variables  are  illustrated  in  Figure  2. 


Equation  (3)  was  rearranged  to  give  an  implicit  equation  for  the 
stretching  coefficient. 

CV  -  h  ln(l+^p  (eC?-D)  (0 

Equation  (4)  was  iteratively  solved  to  determine  the  stretching 
coefficient  in  the  y  direction.  An  algorithm  which  combfned  a  linear 
interpolation,  an  inverse  quadratic  interpolation,  and  a  bisection  was 
used  to  accelerate  the  iterative  scheme  (Ref  21).  The  minimum  grid 
spacing  used  in  equation  (4)  was  the  thickness  of  the  plate. 


A  similar  minimum  grid  spacing  was  used  along  the  complete  length 
of  the  plate  in  the  x  direction.  The  respective  forms  of  equations  (3) 
and  (4)  were  solved  to  obtain  the  stretching  coefficient  in  the  x 
direction.  This  coefficient  was  used  to  exponentially  stretch  the  grid 
from  the  minimum  grid  spacing  at  the  ends  of  the  plate  to  approximately 
one  chord  at  the  outer  boundary.  The  finite  difference  grid  in  the 


physical  plane  is  shown  in  Figure  3  while  an  expandec  view  in  Figure  4 
shows  the  crid  near  the  plate.  The  correspondi ng  computational  domain 
is  shown  in  Figure  5. 

B.  Equations  of  Motion 

The  Navier-Stokes  equations  could  have  been  formulated  in  several 
reference  systems.  The  reference  system  was  selected  to  simplify  the 
overall  approach  to  the  problem.  Two  reference  systems  looked  appealing 
for  application  to  this  problem.  The  first  was  a  rotating  reference 
system  attached  to  the  plate.  The  second  formulation  used  a  mixed 
reference  system  based  on  an  inertial  system. 

The  flow  field  was  calculated  using  the  two  dimensional  time 
dependent  incompressible  Navier-Stokes  equations  in  primitive  variables 
(Refs  22  and  23).  The  problem  was  formulatelJn  a  body  fixed  coordinate 
system  which  has  proven  to  be  a  powerful  tool  in  applying  numerical 
methods  to  practical  boundary  value  problems  involving  partial 
differential  equations  (Ref  6).  The  surface  of  the  plate  was 
represented  as  coordinate  lines  which  reduced  the  difficulties 
associated  with  numerically  specifying  boundary  conditions  in  finite 
difference  methods.  Formulating  the  mathematical  problem  in  a 
coordinate  system  fixed  to  a  body  undergoing  time  dependent  motion 
conveniently  eliminated  the  time  dependent  grid  transformation  terms. 

The  velocities  used  in  the  Navier-Stokes  equations  can  be  cast  in 
several  reference  systems  (Ref  24).  In  most  formulations  the  velocities 
are  cast  in  the  rotating  system  which  means  the  velocities  appear  as  if 
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one  is  movir.c  v.'ith  the  body  (Ref  25).  Hegr.e  (Ref  26)  usee  this  approach 
in  his  pitching  airfoil  problem.  Alternately,  the  velocities  can  be 
cast  in  an  inertial  reference  system  while  the  coordinates  are  ir.  the 
rotating  system.  This  formulation  is  celled  a  "mixed"  system  in  this 
study. 

The  two  velocity  formulations  were  equally  appealing  at  the 
beginning  of  this  study.  A  rotating  reference  system  fixed  to  the  plate 
was  found  to  be  unsatisfactory  because  the  solution  diverged  at  fast 
rotation  rates  as  outlined  in  Appendix  A.  The  other  classical  reference 
system  is  a  purely  inertial  system  in  which  the  plate  moves  within  the 
computational  domain.  The  inertial  system  presents  several 
difficulties.  First,  the  boundary  conditions  imposed  on  the  plate  no 
longer  coincide  with  the  location  of  grid  points.  This  increases  the 
required  computing  time  and  decreases  the  accuracy  because  the  boundary 
conditions  must  be  interpolated.  Second,  the  grid  points  are  no  longer 
properly  clustered  to  resolve  the  large  velocity  gradients  near  the 
plate.  Finally,  with  enough  time,  the  plate  will  exit  the  computational 
doma i n . 

The  Navier-Stokes  equations  were  formulated  in  a  "mixed"  system. 
The  spatial  coordinates  were  formulated  in  a  reference  system  attached 
to  the  plate  (Ref  27).  The  "mixed"  system  retained  the  advantage  of  the 
body  fixed  system  of  calculating  the  grid  only  one  time  and  having  it 
properly  stretched  throughout  the  calculation.  The  flow  velocities  were 
calculated  in  an  inertial  system  to  avoid  the  large  velocities 
introduced  by  the  body  fixed  system.  There  was  no  flow  near  the  outer 
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boundary  in  the  "mixed"  system.  The  velocities  used  in  the  "rixec" 
system,  called  contravariant  velocities  (Ref  21),  were  the  components  of 
the  inertial  velocities  in  the  direction  of  the  coordinate  reference 
system  fixed  to  the  plate.  This  system  avoided  the  difficulty  of  the 
plate  exiting  the  computational  domain  usually  encountered  ir,  an 
inertial  reference  system. 


The  time  dependent,  two  dimensional,  incompressible  Navier-Stokes 
equations  cast  in  primitive  variables  in  an  inertial  reference  system 
are  given  by  (Ref  29): 
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The  body  gravity  forces  were  neglected  even  though  the  aerodynamic 
forces  involved  in  this  problem  were  very  small.  The  classical 
arguments  for  neglecting  the  body  gravity  forces  were  still  valid  (Ref 
30).  The  changes  in  pressure  acting  on  the  plate  caused  by  differences 
in  height  were  insignificant  compared  to  the  other  forces  acting  in  the 
problem. 


The  following  nondimensional  variables  are  introduced: 


(8) 


15 


V  =  * 

"  c 


u 

u  =  - 


V 

V  =  — 

U 


=  P~P°° 


The  nondimensional  variables  were  substituted  into  equations  (5)  through 
(7)  resulting  in  the  following  set  of  equations: 
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The  nondimensional  symbols  are  understood  and  omitted  for  simplicity. 


The  equations  are  transformed  from  an  inertial  (unprimed)  to  a 
rotating  coordinate  system  (primed)  using  the  following  transformations: 


Transformed  Coordinates 

=  _i’  cos  0(t')  -  j_'  sin  0(t’) 

■  i_’  sin  0(t')  +  j_'  cos  Q(t') 

x  ■  xo(t’)  +  x'  cos  0(t’)  -  y'  sin  0(t') 

y  *  y0(t’)  +  x'  sin  0(t')  +  y'  cos  QCt') 


Inverse  Transformed  Coordinates 

i’=  ±  cos  0(t)  +  sin  0(t)  (-*' 

l'=  -  i  sin  e(t)  +  cos  0(t)  (22) 

x’  =  (x-x0(t))  cos  0(t)  +  (y-yQ(t))  sin  e(t)  (22) 

*  -  (x-xo(t))  sin  0(t)  +  (y-yo(t))  cos  £(t)  (24) 


The  transformed  coordinates  and  inverse  transformed  coordinates  in  the 
inertial  axis  and  rotating  body  axis  are  shown  in  Figure  6. 


The  chain  rule  was  used  to  transform  equations  (5)  through  (7)  from 
the  purely  inertial  reference  system  into  the  rotating  body  inertial 
axis  system  given  below: 

ctut'  +  xtux’  +  ytUy'  +  U^txUt'  +  *xux’  +  -y^y'^ 

+  v(ty  ut,+  XyUx,+  vyuy , )  *=  -  (t^p„  +  x^.px  +  >xPy,) 
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The 

transformation  derivatives,  calculated  ■‘‘ror. 

equations  (*7) 

through 

(24) ,  are: 
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The  dot  over  a  term  denotes  differentiation  with  respect  to  time. 


The  transformation  metrics  were  substituted  into  equations  (25) 
through  (27)  end  the  equations  simplified  yielding: 

u'  +  [  (v'e  -  x  cos  G  -  v  sin  £)  +  u  cos  G  +  v  sin  £  ]u  , 
t  ‘  o  o  x 

+  f(-x'G  +  >:  sin  G  -  y  cos  £)  -  u  sin  £  +  v  cos  0]u  , 

0  0  V 
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rx  y  Re  x  x  y  y 


v'  +  [(v'0  -  x  cos  G  -  sin  £)  +  u  cos  ©  +  v  sin  G]v  , 
too  x 

+  [(-x'G  +  x  sin  0  -  y  cos  0)  -  u  sin  0  +  v  cos  ©]v  , 
oo  y 

-  -  (sin  ©  p  ,  +  cos  0  p  . )  +  (v  ,  ,  +  v  ,  ,) 
m  *x  y  Re  x  x  y  y 


cos  S’  u  ,  -  sin  G  u  ,  +  sin  ©  v  ,  +  cos  0  v  .  =  0 
x  y  x  y 


The  contravariant  velocity  components  defined  as  follows: 


u  *  u  cps  0  +  v  sin  ? 
v  *  *u  sin  6  +  v  cos  6 

were  substituted  into  equations  (37)  through  (3S)  giving: 
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B  =  -0u  -  (y’e  -  (xq  cos  0  +  y  sin  0))  v%. , 

.  *  •  • 

-  (-x  0  -(xq  sin  0  +  yQ  cos  ©))  v  ’ 

The  tildes  and  primes  are  omitted  after  this  point  for  simplicity. 


(45) 


Terms  A  and  B  are  similar  to  the  source  terms  which  are 
destabilizing  in  the  rotating  formulation  as  described  in  Appendix  A. 
However,  all  but  one  of  the  individual  terms  from  A  and  B  will  be  moved 
to  the  left  hand  side  of  the  equations.  There  they  will  be  added  to  the 
coefficients  of  the  velocity  derivatives  and  will  be  upwind  differenced 
improving  the  convergence  of  the  implicit  iteration  scheme. 


Calculating  pressure  in  the  primitive  variable  formulation  was 
difficult  because  pressure  did  not  appear  in  the  continuity  equation  nor 
was  there  a  time  derivative  of  pressure  in  the  governing  equations.  Two 
approaches  have  been  extensively  used  to  compute  the  pressure  which 
employed  the  velocity  divergence,  v._V=D',  as  a  correction  factor.  The 
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■first  methoc,  introduced  by  Chorin  (Ref  31),  ccrputed  pressure  using  er. 
iterative  procedure  with  D'  as  the  correcting  term.  The  expression  is 


(s+1)  (s) 

p  =  p 


4>D' 


('r. 


where  s  denotes  the  iteration  number.  The  second  method  used  a  Poisson 
equation  for  pressure  derived  from  the  divergence  of  the  momentum 
equations  expressed  in  equations  (25)  and  (26)  (Refs  16,  32,  and  33). 

The  simplified  incompressible  Poisson  pressure  equation  is 

.  •>  -> 
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(47) 


This  equation  contains  both  spatial  and  time  derivatives  of  the  velocity 
divergence  as  well  as  the  fluid  velocities  and  pressure.  Further 
simplifications  occurred  when  small  variations  in  the  velocity 
divergence  are  assumed.  In  this  case,  the_-spatial  derivatives  of  D' 
were  set  to  zero  in  equation  (47)  and  the  result  becomes 
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(46) 


This  equation  is  identical  to  the  one  used  by  Hodge  (Ref  18)  for  his 
nonrotationg  calculations.  The  Poisson  pressure  equation  technique  was 
used  to  calculate  the  interior  flow  field  because  of  the  increased 
coupling  which  occurred  in  the  finite  difference  representation.  The 
pressure  on  the  plate  was  calculated  using  the  iterative  technique  of 
equation  (46).  Conservation  of  mass  was  imposed  directly  and  the 
difficulty  of  evaluating  the  Laplacian  of  pressure  on  the  plate  was 
avoided. 
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The  following  expressions  taken  fron;  equations  '1:7),  ( 4£ ) ,  (47), 
anc  (48)  became  the  system  used  for  the  sclutior  cf  i  plate  rotating  in 
an  incompressible  viscous  two  dimensional  flow: 

•  *  • 
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The  governing  equations  were  transformed  to  a  computational  plane. 
The  corresponding  governing  equations,  transformed  using  the 
relationships  described  by  Hodge,  become: 
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At  this  point  the  equations  were  valid  for  an  arbitrary 
transformation.  In  the  exponential  transformation  described  in  Section 
III -A ,  each  spatial  coordinate  in  the  computational  plane  was  a  function 
of  only  one  spatial  variable  in  the  physical  plane.  Therefore,  and 
y^  were  both  set  to  zero  which  significantly  simplified  the  governing 
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A  backwarci-time-central-space  finite  difference  method  similar  tc 
the  one  developed  by  Hodge  (Ref  16)  was  used  tc  solve  equations  (53) 
through  (56).  The  implicit  method  consisted  of  a  linearization  of  the 
equations  and  the  application  of  finite  differences  to  these  equations 
(Refs  34  and  35). 

Standard  finite  difference  expressions  were  used  to  approximate 
equations  (53)  through  (56)  (Ref  36).  First  order  accurate  backward 
differences  were  used  for  the  time  derivatives.  The  first  derivatives 
of  pressure  were  approximated  with  second  order  accurate  central 
differences.  Second  order  central  difference  expressions  were  used  for 
the  second  derivatives  of  velocity  in  the  viscous  terms.  Second  order 
accurate  upwind  differencing  was  used  for  the.-velocity  first  derivatives 
in  the  convective  terms.  Three  point,  one  sided,  backward  or  forward 
finite  differences  were  used  for  the  velocity  first  derivatives  in  the 
viscous  terms.  First  order  upwind  differences  were  usea  next  to  the 
plate  and  the  outer  boundary.  The  grid  spacing  near  the  plate  was 
small,  as  were  the  velocity  gradients  near  the  outer  boundary.  In  both 
cases  the  artificial  viscosity  induced  by  first  order  upwind 
differencing  had  a  negligible  effect.  The  transformation  derivatives 
were  evaluated  analytically  or  with  second  order  central  differences. 
The  velocity  derivatives  in  the  velocity  divergence  forms  were  evaluated 
using  central  differences. 

Difference  equations  were  written  in  the  following  convention  used 


b>  hfccne  (Ref  3t).  Terms  which  contained  unknowns  at  location  (i,j)  art 
fully  subscripted.  The  subscripts,  i  and  j,  anc  superscript,  n,  are 
assumed  when  omitted.  Equation  53,  the  x  momentum  eouation,  becomes: 
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where  the  direction  of  the  upwind  differencing  is  determined  by  the  sign 
of  the  coefficient  multiplying  the  term.  This  equation  is  written  to 
provide  an  expression  for  u... 
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In  an  analogous  manner,  the  transformed  y  momentum  equation  is  xr 
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The  static  pressure  throughout  the  flov;  field  was  calculated  using 
the  Poisson  pressure  equation  (Refs  38  and  39).  Second  order  accurate 
central  differences  were  used  to  approximate  the  second  derivatives  of 
pressure  in  the  transformed  Laplacian  term.  The  firs;,  derivative  of 
pressure  in  the  Laplacian  was  approximated  by  second  order  accurate 
upwind  differences  based  on  the  sign  of  o  and  t.  The  remaining  terms 
were  differenced  in  the  same  manner  as  the  momentum  equations.  The 
requirement  that  the  velocity  divergence  at  the  nth  time  step  be  zero 
was  incorporated  into  the  difference  equation  by  setting  D,n  to  zero. 
The  numerical  non-zero  quantity  D,n”^  was  retained  as  a  correction 
factor.  Equation  (55)  was  written  in  finite  difference  form. 
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Equations  (63)  through  (65)  formed  a  set  of  difference  equations 
for  the  flow  variables  u,  v,  and  p  at  location  (i,j)  in  terms  of  the 
transformation  derivatives  and  the  values  of  the  variables  at 
neighboring  points.  The  linearization  and  aifferencing  methods 
decoupled  the  equations  with  respect  to  the  higher  order  terms  for  the 
unknowns  at  location  (i,j). 


The  large  system  of  finite  difference  equations  was  solved  using  a 
successive-over-relaxation  (SOR)  iteration  technique  given  by 
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where  s  is  the  iterate  number  and  *  designates  the  current  value 
obtained  from  the  difference  equation  (Ref  40).  The  acceleration 
parameters  were  always  1.0  for  the  first  iteration. 


The  pressure  on  the  plate  was  calculated  separately  using  the 
Chorin  iteration  technique  which  avoided  the  difficulty  of  formulating 
one  sided  differences  at  the  plate  for  the  Laplacian  of  pressure. 
Second  order  accurate  upwind  differences  were  used  for  the 
transformation  and  velocity  first  derivatives  in  equation  (56).  The 
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dcceleretior  parameter,  < ,  which  if  related  tc  SGR  iteration  of  the 
Rcisscn  pressure  equation  was  given  by  Hodge  (Ref  15)  as 


(69) 


The  iteration  procedure  followed  a  prescribed  order  at  each  point. 
Only  the  surface  pressure  was  computed  on  body  points.  At  field  points, 
the  u  velocity  component  was  computed  first,  then  the  v  velocity 
component,  and  finally  the  field  pressure  p. 


The  SOR  sweep  directions  were  developed  to  minimize  favoring  any 
one  direction.  Hodge  (Ref  41)  found  the  flow  fields  on  the  upper  and 
lower  surfaces  were  slightly  different  on  a  symmetrical  airfoil  at  zero 
angle  of  attack.  This  was  caused  by  always  sweeping  in  the  same 
direction.  The  sweep  method  used  in  this  sturdy  is  described  in  Appendix 
B. 


C.  Initial  and  Boundary  Conditions 

The  initial  conditions  for  the  Navier-Stokes  equations  solver  were 
modified  depending  on  the  application  of  the  solver.  Simple  initial 
conditions  were  used  when  calculating  the  aerodynamics  of  a  plate  at  a 
fixed  angle  of  attack  (Ref  42).  More  refined  initial  conditions  were 
required  when  the  plate  was  allowed  to  rotate  or  free  fall. 

Zero  velocity  was  used  as  the  initial  condition  for  a  plate  at  a 
fixed  angle  of  attack.  An  impulsive  start  was  modeled  by 
instantaneously  imposing  the  freestream  velocity  on  the  plate.  It  was 
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net  neefcsifry  tc  start  with  an  inviscid  solut’oi  because  the  relatively 
',cv  Reynolds  number  did  rict  present  siy  numerical  stability  problems. 
The  numerical  results  did  not  represent  the  physical  phenomenon  durirg 
the  impulsive  start  because  continuity  was  initially  not  satisfied 
throughout  the  flow  field.  As  the  calculation  evolved,  continuity  was 
eventually  satisfied  throughout  the  flow  field  and  the  numerical 
solution  represented  the  physical  phenomenon. 

When  the  Navier-Stokes  equations  solver  was  used  to  calculate  the 
flow  field  around  a  plate  undergoing  forced  rotation,  the  solution  for  a 
plate  at  a  fixed  angle  of  attack  was  used  as  an  initial  condition.  The 
plate  was  then  forced  to  rotate.  As  long  as  the  plate  was  at  a  small 
angle  of  attack  and  was  not  shedding  vortices,  the  start  up  accurately 
modeled  the  physical  phenomenon. 

The  flow  field  around  a  plate  moving  at  a  fixed  angle  of  attack  was 
used  as  an  initial  condition  for  a  plate  about  to  undergo  fixed  or  free 
fall  autcrotation.  This  represented  the  case  of  a  plate  translating  at 
a  constant  velocity  without  rotation.  The  aerodynamic  forces  produced 
by  the  plate  were  used  as  forcing  functions  in  the  flight  mechanics 
equations.  The  flight  mechanics  equations  were  solved  to  initiate 
autorotation  which  accurately  modeled  a  physical  start  up. 

The  Navier-Stokes  boundary  conditions  used  in  this  research  effort 
were  the  reverse  of  those  used  the  more  classical  body  fixed 
formulation.  When  using  a  body  fixed  formulation,  the  boundary 
conditions  were  usually  freestream  velocity  and  pressure  on  the  outer 


boundary  and  no  flow  on  the  inner  boundary  or  plate.  For  the  inertial 
formulation  used  in  this  study,  the  different  boundary  conditions  were 
equivalent  in  their  respective  reference  systems. 


The  motion  of  the  plate  was  imposed  on  the  inner  boundary  or  plate. 
The  translation  and  rotation  of  the  plate  were  modeled  by  changing  the 
boundary  conditions  on  the  inner  boundary.  The  velocity  imposed  at  each 
grid  point  on  the  plate  was  calculated  using  the  following  equations: 


u,  =  V  (cos  £  ccs  y 
b  eg 


sin  0  sin  y) 


-  0v 


(70) 


v,  =  V  (sin  0  cos  y  -  cos  £  sin  y)  +  ©x  (71) 
b  eg 

The  translational  velocity,  angular  velocity,  glide  path  angle,  and 
rotational  angle  were  all  known  from  the  previous  solution  of  the  flight 
mechanics  equations.  Examples  of  various  inner  boundary  conditions  for 

one  side  of  the  plate  are  shown  in  Figure  7. _ Figure  7a  shows  the  inner 

boundary  conditions  imposed  on  a  plate  moving  down  and  to  the  right  at  a 
it/4  radians  angle  in  the  computational  domain.  This  represents  any 
combination  of  flight  path  angle  and  rotation  angle  which  produces  an 
angle  of  attack  of  tt/4  radians  where  the  angle  of  attack  is  defined  as: 


a  =  G-y 


(72) 


Figure  7b  shows  the  inner  boundary  conditions  for  a  plate  rotating 
without  translation.  The  magnitude  of  the  velocity  in  this  case  is 
purely  a  function  of  the  angular  rotation  of  the  plate.  The  translation 
from  Figure  7a  was  combined  with  the  rotation  in  Figure  7b  to  produce 
the  boundary  conditions  shown  in  Figure  7c. 
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The  velocity  boundary  conditions  or  the  outer  boundary  were  zero 
velocity  ir  either  the  x  or  y  directions.  Trie  pressure  on  the  cuter 
boundary  was  set  equal  tc  the  ambient  pressure  or  c  ncndimensional 
pressure  of  zero.  These  boundary  conditions  closely  correspondec  tc  the 
physical  conditions  in  an  inertial  reference  system.  The  outer  boundary 
conditions  were  not  changed  and  remainec  constant  throughout  the 
calculations.  Neglecting  the  momentum  deficit  had  a  small  effect  on  the 
flow  field  near  the  plate  as  outlined  in  Appendix  A. 

0.  Determination  of  Aerodynamic  Coefficients 

The  flow  field  around  the  plate  was  calculated  by  the  Navier-Stokes 
equations  solver.  The  viscous  and  pressure  forces  imposed  upon  the 
plate  by  this  flow  field  were  summed  to  give  the  aerodynamic  forces 
generated  by  the  plate. 

The  aerodynamic  coefficients  were  calculated  at  each  time  step, 
whereas  for  a  steady  state  solution  they  are  only  calculated  for  the 

final  state  (Ref  43).  Because  the  coefficients  were  calculated  so  often 

in  this  study,  the  method  of  calculation  was  customized  to  this 
particular  geometry  to  reduce  the  computing  requirements. 

The  total  pressure  force  acting  on  the  plate  was  calculated  by 
summing  the  pressure  acting  on  all  the  descrete  segments  of  the  plate. 
The  translational  and  rotational  velocities  of  the  plate  were  specified 
as  boundary  conditions.  The  pressure  distribution  on  the  plate  was 

calculated  using  equation  (65).  The  pressure  differences  around  the 
plate  were  summed  to  give  the  normal  and  axial  forces.  The 


corresponding  moment  was  calculated  b\  multiplying  the  pressure 
differences  by  their  respective  distances  from  the  center  of  gravity. 


The  viscous  forces  acting  on  the  plate  were  calculated  using  the 
following  equation  for  shear  stress: 


where  un  is  the  derivative  of  the  velocity  normal  to  the  surface  of  the 
plate  (Ref  44).  The  normal  velocity  derivative  was  approximated  by  a 
first  order  difference.  The  viscous  normal  and  axial  forces  were 
calculated  by  summing  the  shear  stresses  around  the  plate  and  the  moment 
produced  by  the  viscous  forces  was  calculated  by  multiplying  the  shear 
stress  by  its  corresponding  distance  from  the  center  of  gravity. 


The  viscous  and  pressure  forces  were  combined  according  to 


L  *  (F  +  F  )  cos  a  -  (F  +  F  )  sin  a 

yv  yp  xv  xp 

D  =  (F  +  F  )  sin  a  +  (F  +  F  )  cos  a 

w  vc  xv  vp 


M  =F4p(Ax)  +  IAt  (Ax)  y.. 


giving  the  overall  aerodynamic  forces  produced  by  the  plate.  The 
aerodynamic  forces  were  produced  primarily  by  the  pressure  differences 
around  the  plate. 


31 


SECTION  IV 


FLIGHT  MECHANICS  SOLVER 


A  three  degree  of  freedom  flight  mechanics  equations  solver  was 
developed  to  predict  the  flight  path  of  a  free  felling,  autorotating, 
flat  plate.  The  three  flight  mechanics  equations  were  rewritten  as  four 
coupled  first  order  ordinary  differential  equations.  The  equations  were 
solved  using  an  open  Adams  formula  to  obtain  the  velocity  of  the  center 
of  gravity  of  the  plate,  the  flight  path  angle,  the  rate  of  rotation, 
and  the  rotational  angle. 

The  general  unsteady  motion  of  a  body  in  three  dimensional  space  is 
described  by  a  vector  force  and  moment  equation.  This  produces  six 
scalar  ordinary  differential  equations.  Much  of  the  complexity  of  these 
equations  results  from  the  inclusion  of  the__rotation  of  the  earth. and 
its  curvature  in  the  mathematical  model.  Assuming  the  earth  is  flat, 
the  absence  of  gyroscopic  effects  and  planar  motion  reduces  the  system 
of  equations  to  three  coupled  differential  equations  which  model  the 
motion  of  a  free  falling,  autorotating  plate  in  a  plane  (Ref  45). 

Two  alternative  reference  systems  are  available  for  the  dynamical 
force  equations:  the  wind  axes  and  the  body  axes.  Both  are  used  in 
practice  with  no  overriding  advantages  to  either  system  (Ref  46).  The 
equations  used  in  this  research  effort  were  cast  in  the  wind  axes 
reference  system. 

Experimental  data  and  previous  analyses  indicated  the  motion  of  the 
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autc rota ting  plate  could  be  modeled  accurately  using  the  following  three 
degrees  of  freedon,; 

1)  Horizontal  translation,  x 

2)  Vertical  translation,  y 

3)  Rotation  about  the  center  of  gravity,  e 

The  reference  system  for  these  three  degrees  of  freedom  is  shown  in 
Figure  8  (Ref  47).  The  positive  x  direction  was  to  the  right  while  the 
positive  y  direction  was  up.  The  angle  of  rotation,  like  all  angles  in 
this  problem,  is  positive  in  a  counter-clockwise  direction.  This 
conforms  to  a  standard  right  hand  rule  reference  system. 

The  equations  of  motion  modeling  these  three  degrees  of  freedom 

are: 

-mV  =  D  +  W  sin  Y  (77) 

nVY  =  L-W  cos  v  —  ( 78) 

Iv  £'  =  M  (79) 

This  was  an  initial  value  problem  involving  a  system  of  three  nonlinear 
ordinary  differential  equations  (two  first  order  and  one  second  order). 
Several  methods  exist  for  solving  systems  of  first  order  equations, 
therefore,  the  second  order  equation  was  reduced  to  two  first  order 
equations. 

I  q  *  K  (80) 

y 

®  “  <3  (81) 

Equation  (79)  was  reformulated  using  equations  (80)  and  (81)  to  give 
equations  (82).  Equations  (77)  and  (78)  were  rearranged  to  give 
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equations  (83)  end  (84). 


c-  =  f~  '82- 

y 

v  =  -CIH-W  sin  y)  (£3 ) 

m 

=  Q--y  cos  y}  (84) 

mV 

Several  implicit,  or  closed,  and  several  explicit,  or  open, 
methods  are  available  for  solving  systems  of  equations  (Ref  48).  The 
closed  methods  offer  improved  accuracy  over  the  open  methods,  however, 
they  require  multiple  solutions  of  the  Navier-Stokes  equations  at  each 
time  step.  During  the  formulation  of  this  problem,  it  was  estimated  the 
Navier-Stokes  equation  solver  would  consume  at  least  90%  of  the 
computing  time  while  the  flight  mechanics  equations  solver  would  consume 
less  than  10%.  Based  on  this  assumption,  an  open  initial  value  problem 
solver  was  selected  as  the  most  computationally  efficient  method.  A 
predictor-corrector  scheme  was  ruled  out  because  it  required  two 
solutions  of  the  Navier-Stokes  equations  at  each  time  step,  thus  almost 
doubling  the  computing  time. 

The  simplest  explicit  method  is  Euler's  method.  This  method 
requires  one  initial  condition  for  each  variable  and  one  evaluation  of 
the  function  at  each  time  step  to  advance  to  the  next  time  step.  The 
following  set  of  equations  represents  Euler's  method: 

Vn+1  =  Vn  +  AtV 
’fn+1  *  yn  +  Aty 
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(85) 

(86) 


Euler's  method  is  first  order  accurate  in  time  which  is  consistent  with 
the  time  differencing  used  in  the  Navier-Stokes  equations  solver. 

Few  additional  computer  resources  are  required  to  increase  the 
accuracy  of  the  method  by  using  higher  order  Adams-Bashforth  formulae 
(Ref  48).  These  methods  are  described  in  Appendix  D.  Equations  (81) 
through  (84)  were  nondimensionalized  using  the  same  operators  that  were 
used  in  the  Navier-Stokes  equations  solver.  The  equations  are  cast  in 
the  Adams-Bashforth  formulae 

n+l  =  ,.n  C£_  r  aCOEF.  (D.  +  V  sin  Y  )  • 

“  WV  1=1  1  1  1 

CD 

,  n  L.  -  V  cos  y 

yn+i  =  Yn  +  At  £S_  .1  aCOEF.  (— - ^ - )  (90) 

co  -  i 

qn+1  =  q11  +  At  iI1  ACOEFi  (Q1) 

00 

n 

en+1  =  en  +  At  iI  ACOEF,  q  (C9) 

where  n  represents  the  order  of  the  solution  technique.  A  list  of  the 
coefficients  used  in  the  summation  are  included  in  Appendix  D.  If  the 
order  of  the  equations  is  set  to  one,  equations  (85)  through  (88)  reduce 
to  equations  (89)  through  (92).  Increasing  the  order  of  the  flight 
mechanics  equations  solver  had  an  insignificant  effect  on  the  overall 
storage  requirements  and  computational  times.  An  Euler  solution 
required  storing  20  pieces  of  information.  Each  increase  in  the  order 
of  the  solver  required  saving  an  additional  10  pieces  of  information. 


Therefore,  e  sixth  order  method  required  saving  7G  pieces  of  information 
as  opposed  to  20  for  a  first  order  solver.  This  was  insignificant 
compared  to  the  thousands  of  storage  locations  required  for  the 
Navier-Stokes  equations  solver.  A  comparison  between  different  order 
Adams-Bashforth  formulae  was  made  and  is  discussed  in  Section  VI. 

The  initial  conditions  necessary  to  start  the  higher  order  methods 
are  not  known.  A  common  technique  is  to  start  with  Euler's  method  and 
increase  the  order  of  the  method  at  each  time  step  until  the  desired 
order  is  obtained.  This  in  essence  generates  the  initial  conditions 
necessary  to  start  the  higher  order  techniques.  In  this  study,  the  past 
initial  conditions  necessary  to  start  the  higher  order  methods  were 
obtained  by  linear  extrapolation  backward  in  time  from  the  one  set  of 
initial  conditions  which  were  specified. 

Everything  on  the  right  hand  side  of  equations  (£9)  through  (92)  is 
known.  The  velocity,  glide  path  angle,  rotational  angle,  and  rotational 
velocity  were  all  known  from  previous  solutions  of  the  flight  mechanics 
equations.  The  aerodynamic  forces  (lift,  drag  and  moment)  are  known 
from  the  present  and  previous  time  step  solutions  of  the  Navier-Stokes 
equations.  The  aerodynamic  forces  and  plate  characteristics  were  all 
referenced  to  a  plate  with  a  span  of  1.0. 
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COUPLED  SYSTEMS 


The  Navier-Stokes  equations  solver  and  the  flight  mechanics 
equations  solver  were  solved  simultaneously.  Initial  conditions  were 
specified  for  both  sets  of  equations.  The  aerodynamic  forces  produced 
by  the  plate  were  calculated  by  the  Navier-Stokes  equations  solver  and 
used  as  forcing  functions  in  the  flight  mechanics  equations  solver.  The 
flight  mechanics  equations  were  solved  to  obtain  the  translation  and 
rotation  of  the  plate  for  the  next  time  step.  The  updated  translation 
and  rotation  were  used  as  new  boundary  conditions  in  the  Navier-Stokes 
equations  solver.  The  iterative  procedure  was  repeated  thousands  of 
times,  therefore,  a  solution  dependent,  variable  time  step  was 

introduced  into  the  computational  procedure  to  improve  the  computational 
efficiency. 

Zero  velocity  in  the  computational  domain  was  used  as  the  initial 
condition  in  the  Navier-Stokes  equations  solver.  The  plate  was 

impulsively  accelerated  to  some  prespecified  velocity  at  a  fixed  angle 
of  attack.  A  relatively  large  time  step  was  used  to  converge  to  the 
steady  state  solution  with  the  minimum  expenditure  of  computer  time. 
Therefore  the  flow  field  which  evolved  during  start  up  did  not 
accurately  represent  physics.  The  solution  at  a  fixed  angle  of  attack 
included  periodic  shedding  of  vortices  when  the  angle  of  attack  was  more 
than  a  few  degrees.  The  intent  at  this  point  was  not  to  accurately 
model  the  impulsive  start  up,  but  rather  to  obtain  a  fully  developed 
flow  field  which  could  be  used  to  start  the  rotating  calculations.  The 


aerodynamic  forces  were  calculated  frcr  the  pre ssi.>'e  anc  velocity  field 
arour.d  the  plate. 

The  initial  conoitions  for  the  flight  mechanics  equations  were 
specified  at  the  beginning  of  the  calculation.  The  initial  conditions 
for  previous  time  steps  necessary  to  start  the  higher  order  solvers  were 
linearly  extrapolated  backwards  in  time.  The  aerodynamic  forces 
produced  by  the  plate  were  entered  into  the  flight  mechanics  i=s  boundary 
conditions  in  equations  (89)  through  (92).  The  same  time  step  was  used 
in  the  flight  mechanics  equations  solver  as  was  used  in  the 
Navier-Stokes  equations  solver.  The  flight  mechanics  equations  were 
solved  for  the  next  time  step.  This  solution  produced  a  new  velocity, 
glide  path  angle,  rotation  angle,  and  rotational  velocity. 

These  conditions  were  used  in  equations  (71)  and  (72)  to  obtain  new 
inner  boundary  conditions  for  the  Navier-Stokes  equations  solver.  The 
outer  boundary  conditions  of  no  flow  and  ambient  pressure  remained 
unchanged.  A  new  flow  field  and  the  corresponding  aerodynamic 
coefficients  were  calculated.  At  this  point  the  solution  procedure  was 
started.  The  motion  of  the  free  falling  plate  was  calculated  by 
repeating  the  process  just  outlined. 

The  time  steps  required  to  obtain  an  accurate  Navier-Stokes 
solution  and  an  accurate  flight  mechanics  solution  were  not.  the  same. 
The  Navier-Stokes  equations  solver  was  an  implicit  scheme.  A  stability 
analysis  of  an  implicit  technique  used  to  solve  a  model  linear  partial 
differential  equation  with  the  same  character  as  the  Navier-Stokes 
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ecu  tions  showed  the  numerical  method  tc  be  unconditionally  stable  (Re: 
36; .  However,  when  working  with  the  complete  system  of  Navier-Stokes 
ecueticns,  Hegna  (Ref  37)  found  the  existence  or  an  upper  limit  on  the 
time  step.  The  flight  mechanics  equations  were  stable  for  any  time 
step.  Therefore,  the  time  step  was  selected  to  meet  the  limitations  of 
the  Navier-Stokes  equations  solver. 

The  Navier-Stokes  equations  solver  was  found  to  tolerate  only  a 
certain  unquantifiable  level  of  change  in  the  boundary  conditions  from 
time  step  to  time  step.  The  solver  still  met  the  error  criteria 
specified  for  convergence  within  an  acceptable  number  of  iterations. 
This  occured  because  the  Navier-Stokes  equations  were  not  solved 
simultaneously .  The  updated  velocity  boundary  conditions  on  the  plate 
introduced  large  changes  in  u  and/or  v  near  the  body.  These  changes  in 
velocity  caused  large  changes  in  pressure  on  the  plate  through  the 
velocity  divergence  term  in  equation  (46).  The  large  pressure  changes 
in  turn  influenced  v  at  the  points  adjacent  to  the  plate  because  cf  the 
coupling.  This  interaction  between  equations  made  convergence 

difficult  for  large  changes  in  the  boundary  conditions  on  the  plate. 
Numerical  experimentation  was  used  to  establish  the  tolerable  level  of 
change. 

The  maximum  allowable  change  in  the  boundary  conditions  was  usee  as 
the  basis  for  c  solution  dependent  time  step.  The  maximum  change  in  the 
boundary  condition  was  asssumed  to  be  a  result  of  the  rate  of  rotation 
of  the  plate.  A  heuristic  method  was  developed  around  the  following 
equation: 
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The  constant  for  a  variable  time  step,  C„,  was  established  through 
numerical  experimentation  to  a  value  which  allowed  the  SOR  procedure  to 
meet  the  specified  convergence  criteria.  This  allowed  approximately  the 
same  amount  of  rotation  during  each  time  step  regardless  of  the  rate  of 
rotation.  When  the  plate  was  rotating  rapidly,  a  relatively  small  time 
step  was  used  resulting  in  a  certain  amount  of  change  in  ,,te  boundary 
conditions  between  calculations.  When  the  plate  was  rotating  slowly,  a 
relatively  large  time  step  was  used  resulting  in  approximately  the  same 
amount  of  change  in  the  boundary  conditions  as  the  rapidly  rotating 
case.  Without  the  solution  dependent  time  step,  a  very  small  time  step 
would  have  been  required  to  accurately  calculate  the  most  rapid  rotation 
thereby  resulting  in  inefficiency  at  slower  rotation  rates. 

The  higher  order  Adam-Bashforth  formulae  used  to  solve  the  flight 
mechanics  equations  were  derived  for  constant  time  steps-.  The  first 
order  formula  or  Euler's  method  was  the  only  scheme  which  was  used  with 
the  variable  time  step.  Higher  order  methods  with  variable  time  steps 
exist,  however,  their  derivation  is  rather  complex. 
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section  v: 

DISCUSSION  0-  RESULTS 


Several  major  stepping  stones  were  investigated  cr.  the  way  to 
coupling  the  Navier-Stokes  equations  with  the  flight  mechanics 
equations.  A  building  block  approach  was  used  to  validate  the  overall 
approach  and  the  computer  code.  Four  different  cases  were  investigated 
in  this  study. 

(A)  The  Navier-Stokes  equations  solver  was  used  to  calculate  the 
flow  field  around  a  plate  at  several  static  angles  of  attack. 

(B)  The  effect  of  forcing  the  plate  to  rotate  was  also  investigated 
using  the  Navier-Stokes  equations  solver. 

(C)  The  Navier-Stokes  equations  were  coupled  with  a  one  degree  of 
freedom  equation  of  motion  and  solved  to  model  the  pinned  autorotation 
of  a  plate. 

(D)  Finally  the  flight  path  of  a  free  falling,  autorotating,  flat 
plate  was  predicted  by  a  numerical  method  which  coupled  the  two 
dimensional  incompressible  Navier-Stokes  equations  with  the  three 
degrees  of  freedom  equations  of  motion. 

The  finite  difference  grid  which  modeled  the  flow  field  was 
constructed  to  obtain  approximately  the  same  resolution  in  the  x  and  y 
directions.  The  minimum  and  maximum  grid  spacing  determined  the  total 
number  of  grid  points  required  to  cover  the  flow  field.  Numerical 
experimentation  was  used  to  select  a  grid  with  a  maximum  number  of 
points  compatible  with  calculations  on  the  CYBER  750.  The  minimum  grid 
spacing  in  the  y  direction,  which  defined  the  thickness  of  the  plate. 
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was  chore  or  C . 0156  feet  cr,  a  plate  with  a  C.25  feet  chord.  The 

minimum  grid  spacing  and  the  requirement  for  an  outer  spacing  of 
approximately  one  chord  dictated  a  grid  with  5A  points  in  the  y 
direction.  An  exponential  stretching  coefficient  of  2.988  produced  a 
maximum  grid  spacing  of  1.110  chords  on  the  outer  boundary  in  the  y 
direction.  The  physical  domain  extended  10.03125  chords  in  the  y 
direction,  or  2.508  feet  for  a  plate  with  a  0.25  feet  chord. 

The  same  minimum  grid  spacing  was  used  in  the  x  direction.  Sixteen 
grid  points  were  evenly  spaced  along  the  chord  of  the  plate  at  6.25% 
chord  intervals.  The  grid  was  stretched  from  the  end  of  the  plate  to 
the  outer  boundary  using  an  exponential  coefficient  of  3.157.  This 
required  65  points  to  produce  an  outer  grid  spacing  of  1.287.  The 
physical  domain  extended  10.5  chords  in  the  x  direction  or  2.625  feet 
for  a  plate  with  a  0.25  feet  chord.  Hence  the  complete  grid  was 
comprised  of  3510  points. 

The  finite  difference  grid  was  generated  by  piecing  together 
analytical  functions.  This  presented  the  opportunity  to  use  analytical 
transformation  derivatives  as  well  as  the  more  commonly  used  numerical 
derivatives.  The  flow  field  around  a  plate  at  t/2  radians  angle  of 
attack  was  calculated  using  analytical  and  numerical  derivatives.  The 
analytical  transformation  derivatives  produced  a  drag  coefficient  five 
drag  counts  (0^=0.0005)  higher  than  the  numerical  derivatives. 
Analytical  derivatives  were  used  throughout  the  study  because  they  were 
thought  to  be  more  accurate. 


The  Reynolds  number  cf  a  plate  with  a  0.25  feet  chord  moving  et  a 
freestream  velocity  of  3.0  feet  per  second  on  a  standard  day  at  sea 
level  is  approximately  5000.  The  Navier-Stokes  equations  were  not  mess 
averaged  because  of  the  relatively  small  Reynolds  number  (Ref  49).  No 
turbulence  modeling  was  included  in  the  equations. 

Implicit  methods  were  shown  by  linear  stability  analysis  to  be 
unconditionally  stable  (Ref  50).  Despite  this,  Hegna  (Ref  51)  found 
convergence  of  the  S0R  solver  required  a  nondimensional  time  step  less 
than  0.001  for  flow  around  an  airfoil  at  a  Reynolds  number  of 
approximately  170,000.  The  relatively  low  Reynolds  number  flow  in  this 
study  produced  a  system  of  equations  which  could  be  solved  using  a 
larger  time  step.  Through  numerical  experimentation,  a  time  step  of 
0.01  was  found  to  converge  and  was  used  for  all  calculations  at  fixed 
angles  of  attack. 

There  was  a  tradeoff  in  the  S0R  procedure  between  the  size  of  the 
time  step  and  the  number  of  iterations  performed  at  each  time  step. 
Increasing  the  time  step  and  increasing  the  number  of  iterations  at  each 
time  step  can  have  the  same  effect  as  decreasing  the  time  step  and 
decreasing  the  number  of  iterations.  In  this  study  a  maximum  of  four 
S0R  iterations  were  used  at  each  time  step.  An  even  number  of 
iterations  was  always  used  in  an  attempt  to  avoid  using  one  sweep 
direction  more  than  another. 

The  same  S0R  acceleration  parameters  were  used  throughout  the 
study.  The  parameters  had  values  of  1.0  for  the  velocity  components. 
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1.9  ^or  the  surface  pressure,  and  1.1  fcr  the  field  pressure.  These 

parameters  './ere  determined  by  previous  researchers  and  their  influence 
on  convergence  was  not  investigated  in  this  stuoy  (Refs  16  enc  37). 

The  convergence  criterion  for  the  S0R  iterations  required  the 
maximum  change  in  the  relative  magnitudes  of  u,v,  and  p  for  all 
locations  be  less  than  1%.  This  criterion  occasionally  was  not 

satisfied,  however,  previous  researchers  tolerated  maximum  changes  as 
large  as  5%  (Ref  37).  The  convergence  criterion  was  maintained  at  a 
relatively  low  level  throughout  this  study  because  of  the  strong  time 
dependent  nature  of  the  profem. 

Several  methods  of  SOR  sweeping  were  investigated.  The  method 
described  in  Appendix  B  was  found  to  be  the  best  combination  of 

efficiency  and  accuracy.  This  method,  while  not  favoring  any  direction, 

maintained  good  numerical  efficiency  because  it  did  not  require  ar 
excessive  amount  of  logic  to  execute. 

The  movement  of  the  plate  was  calculated  using  first  order  and 
sixth  order  flight  mechanics  equations  solvers.  There  was  virtually  no 
difference  in  the  solutions  because  the  small  time  step  required  to 
solve  the  Navier-Stokes  equations  minimized  the  error  in  the  first  order 
flight  mechanics  equations  to  an  acceptable  level. 

The  numerical  technique  required  approximately  0.645  seconds  of 
central  processing  time  on  the  CYBER  750  computer  to  advance  a  solution 
one  time  step.  This  equated  to  approximately  2.582  seconds  of  computing 
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time  per  SOT.  iteration.  The  computing  time  for  each  SOT  iteration 
included  approximately  0.001  seconds  to  solve  the  flight  mechanics 
equations.  This  means  98.5  percent  of  the  computing  requirements 
support  the  solution  of  the  Navier-Stokes  equations. 

The  numerical  procedure  was  executed  on  the  CRAY-1  computer  to 
reduce  the  overall  computing  time.  The  computer  code  was  not  modified 
to  incorporate  the  vectorized  programing  techniques  available  on  the 
CRAY-1.  The  computing  time  on  the  CRAY-1  was  slightly  more  than  a  fifth 
of  the  computing  time  required  on  the  CYBER  750  to  perform  identical 
calculations. 

A.  Static  Angle  of  Attack 

The  static  lift  curve  for  a  flat  plate,  shown  in  Figure  9,  was 
calculated  using  the  Navier-Stokes  equations  solver.  A  flow  field  was 
considered  to  achieve  steady  state  by  previous  researchers  when  the 
aerodynamic  coefficients  changed  less  than  1%  in  one  unit  of 
nondimensional  time  (Ref  52).  The  maximum  steady  state  lift  coefficient 
was  approximately  0.44  at  0.15  radians  angle  of  attack.  The  flow 
separated  from  the  upper  surface  of  the  plate  at  a  relatively  small 
angle  of  attack.  Above  0.15  radians  angle  of  attack  the  flow  field  was 
unsteady  and  a  band  of  lift  coefficients  were  plotted  on  the  static  lift 
curve. 

The  instantaneous  streamlines  around  a  plate  at  0.4  radians  angle 
of  attack  can  be  seen  in  Figure  10.  The  streamlines  were  calculated  by 
integrating: 
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using  a  first  order  method.  A  vortex  can  be  seer,  behind  the  plate  in 
Figure  10  (Ref  53).  At  t=5.2  a  vortex  is  just  leaving  the  plate.  It  is 
swept  downstream  until  finally,  at  t=6.0,  another  vortex  is  forming 
above  the  plate.  The  Strouhal  number  for  this  case  was  approximately 
fc/U  =0.34.  The  Strouhal  number,  however,  is  reduced  to  0.13  when  the 

CD 

frontal  area  of  the  plate  was  used  as  the  reference  length.  The 
unsteady  lift  coefficient  produced  by  a  plate  at  0.4  radians  angle  of 
attack  is  illustrated  in  Figure  11. 

The  flow  field  produced  by  a  plate  at  tt/2  radians  angle  of  attack 
was  calculated.  The  symmetric  vortex  formation  behind  the  plate  can  be 
seen  in  the  Figure  12.  Experiment  indicates  that  a  vortex  street  should 
have  been  formed  behind  the  plate  (Refs  54  and  55),  however,  it  was  not 
predicted  by  the  numerical  calculation.  An  instability  of  the  magnitude 
necessary  to  initiate  asymmetric  shedding  was  not  present  in  the 
numerical  method.  No  attempt  was  made  to  initiate  the  asymmetry  because 
it  was  beyond  the  scope  of  the  present  study.  The  drag  coefficient 
produced  by  the  plate  at  tt/2  angle  of  attack  was  1.494.  This  compares 
well  with  the  experimental  data  of  1.4  obtained  for  a  flat  plate  with 
asymmetric  shedding  and  1.6  obtained  for  a  flat  plate  with  a  splitter 
plate  behind  it  (Ref  56). 

6.  Forced  Rotation 

From  this  point  to  the  end  of  this  section,  the  flow  fields 
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presented  combine  translation  and  rotation.  Steady  motions  can  usually 
be  presented  in  a  body  fixed  reference  system  in  which  the  flow  becorr.es 
steady.  In  unsteady  flows  no  preferred  reference  frame  exists, 
therefore,  the  reference  frame  is  selected  on  the  basis  of  other 
criteria  as  discussed  in  Appendix  C.  A  reference  frame  fixed  to  the 
plate  in  translation  in  which  the  plate  rotates  will  be  used  throughout 
this  section.  This  frame  was  selected  because  of  its  analogy  to 
instantaneous  streamlines  of  a  plate  rotating  in  a  wind  tunnel. 

A  solution  dependent  time  step  was  developed  to  improve  the 
computational  efficiency  of  the  rotating  calculations.  The  constant  in 
equation  (92),  which  was  used  to  determine  the  maximum  allowable  change, 
was  derived  from  numerical  experimentation.  The  constant  was  selected 
to  reduce  the  number  of  calculations  required  for  a  solution  while 
maintaining  the  specified  error  criterion  and  using  a  maximum  of  four 
SOR  iterations.  A  maximum  allowable  change  constant  of  Cc=0.005  was 
used  for  most  calculations.  In  addition,  the  time  step  was  never  larger 
than  0.01. 

The  flow  field  around  a  plate  undergoing  forced  rotation  was 
calculated  at  three  nondimensional  rotational  velocities.  The  plate  was 
initially  fixed  in  a  uniform  freestream  flow  at  zero  angle  of  attack. 
It  was  subjected  to  a  constant  angular  acceleration  until  it  reached  the 
desired  rotational  velocity.  First,  the  plate  was  forced  to  rotate 
counterclockwise  at  a  nondimensional  rotational  velocity  (uc/U^)  of  one. 
At  this  rate  of  rotation  the  translational  velocity  everywhere  on  the 
plate  was  greater  than  the  rotational  velocity.  This  meant  vortices 
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formed  at  the  encpoints  cf  the  plate  were  always  shec  downstream  of  the 
plate.  The  instantaneous  streamlines  for  this  case  are  shown  in  Figure 
12. 

The  rate  of  rotation  was  increased  to  an  angular  velocity  of  two. 
The  rotational  velocity  of  the  end  of  the  plate  was  equal  to  the 
translational  velocity  of  the  center  of  gravity  of  the  plate.  This 
meant  that  when  the  plate  was  at  an  angle  of  attack  of  m/2  or  3tt/2  there 
was  no  relative  flow  around  the  retreating  end  of  the  plate  and  thus  no 
vortex  was  shed.  This  can  be  seen  in  Figure  14. 

Finally,  the  rotational  rate  was  increased  to  an  angular  velocity 
of  four.  At  this  rate  of  rotation,  the  linear  velocity  of  the  ends  of 
the  plate  induced  by  the  rotation  was  twice  that  of  the  center  of 
gravity  of  the  plate.  This  meant  a  vortex  was  shed  upstream  of  the 
retreating  end  of  the  plate  as  shown  by  the  instantaneous  streamlines  in 
Figure  15. 

The  aerodynamic  coefficients  produced  by  a  plate  undergoing  forced 
rotation  at  three  different  rates  are  shown  in  Figures  16,  17,  and  18. 
In  these  plots,  angle  of  attack  and  rotation  angle  are  synonymous..  The 
maximum  lift  coefficients  were  approximately  3.2,  5.8,  and  9.5  for 
nondimensional  rotational  rates  of  1.0,  2.0,  and  4.0,  respectively. 
These  values  are  considerably  larger  than  the  maximum  lift  coefficient 
of  approximately  1.2  produced  by  a  flat  plate  at  a  static  angle  of 


attack. 


Rotating  the  plate  caused  a  shift  in  the  phase  '-elatior.ship  among 
the  aerodynamic  coefficients.  The  maximum  drag  coefficient  was  found  to 
leg  the  maximum  lift  coefficient  by  1.1,  1.0,  and  0.75  radians  angle  of 
attack  for  nondimensional  angular  velocities  of  1.0,  2.0,  and  4.0 
respectively.  At  a  rotational  rate  of  one,  the  maximum,  moment 
coefficient  slightly  lagged  the  maximum  lift  coefficient.  At  a 
rotational  rate  of  two  the  phase  lag  increased  to  approximately  0.3 
radians.  The  maximum  moment  lagged  the  maximum  lift  by  approximately 
0.5  radians  at  a  rotational  velocity  of  four.  Visual  integration  of  the 
area  under  the  moment  curve  over  a  revolution  indicated  the  slowest  rate 
of  rotation  produced  an  average  moment  which  would  accelerate  the  rate 
of  rotation.  The  plate  rotating  at  an  angular  velocity  of  two  generated 
an  integrated  moment  which  was  almost  neutral  while  the  most  rapidly 
rotating  case  produced  an  overall  retarding  moment.  This  analysis  of 
the  moment  suggested  the  plate  would  autorotate  at  a  mean  rotational 
velocity  of  approximately  two. 

C.  Pinned  Autorotation 

A  calculation  was  performed  for  a  plate  undergoing  pinned 
autorotation.  This  was  analogous  to  a  plate  being  held  at  its  center  of 
gravity  by  frictionless  bearings  in  a  wind  tunnel.  Pinned  autorotation 
was  numerically  modeled  by  eliminating  equations  (88)  and  (89)  in  the 
flight  mechanics  equations.  The  only  permissible  degree  of  freedom  was 
rotation  which  meant  only  equations  (94)  and  (95)  were  solved  to 
determine  the  motion  of  the  plate. 

The  flow  field  around  a  plate  moving  at  sea  level  in  a  freestream 


velocity  of  3.G  ft/sec  et  e  fixed  ancle  cf  attach  of  0.1  recurs  was 
calculated  and  used  as  ar  initial  condition  fer  rotation.  The  plate 
with  the  physical  characteristics  displayed  in  Table  I  was  "released" 
and  allowed  to  autorotate  by  solving  equations  (94)  end  (95). 


Table  I.  Physical  Characteristics  of  the  Plate 


Chord 

3  inches 

Span 

8  inches 

Aspect  Ratio 

2.67 

Weight 

7.02  x  10'3  lbs 

Moment  of  Inertia 

1.13  x  10‘6  lb* 

The  plate  began  rotating  and  established  autorotation  at  a  mean 
nondimensional  rate  of  rotation  of  1.58.  This  rate  of-  rotation  was 
within  the  range  suggested  by  analyzing  the  moment  in  the  forced 
rotation  studies.  The  aerodynamic  coefficients  produced  by  the  plate 
undergoing  pinned  autorotation  are  shown  in  Figure  19.  In  this  plot, 
angle  of  attack  and  rotation  angle  are  synonymous.  The  magnitude  of  the 
coefficients  are  comparable  to  those  produced  by  the  plate  undergoing 
forced  rotation  at  a  nondimensional  rotational  velocity  of  2.0.  The 
maximum  drag  coefficient  lags  the  maximum  lift  coefficient  by 
approximately  0.70  radians.  The  maximum  moment  coefficient  is  in  phase 
with  the  maximum  lift  coefficient.  The  numerically  integrated  area 
under  the  moment  curve  was  zero  for  a  complete  revolution  after 


50 


ecu'il  ibriun.  an  terete*  ion  was  established. 

The  rate  of  rotation  varied  considerably  during  pinnec  eutorotation 
because  of  the  relatively  small  mass  c*"  the  plate  used  in  this  study. 
E.H.  Smith  (Ref  14)  found  that  the  rate  of  rotation  of  a  relatively 
dense  plate  undergoing  pinned  autorotation  in  a  wind  tunnel  was  almost 
constant.  Lugt  assumed  a  constant  rate  of  rotation  in  his  numerical 
research  (Ref  57).  The  relatively  smell  mass  of  the  plate  in  this  study 
produced  the  variations  in  the  rate  of  rotation  shown  in  Figure  20. 


D.  Free  Fall  Autorotation 

Finally  the  plate  with  the  physical  characteristics  presented  in 
Table  I  was  allowed  to  undergo  free  fall  autorntation.  The  initial 
conditions  for  the  plate  were  moving  down  and  to  the  right  with  a 
velocity  of  3.0  ft/sec  at  a  glide  path  angle  of  -n/4  radians.  The  plate 
was  fixed  at  a  rotation  angle  of  -tr/4+0.1  radians  which  resulted  ir,  an 
angle  of  attack  of  0.1  radians.  The  plate  was  "released"  and  allowed  to 
free  fall  by  solving  the  complete  set  of  flight  mechanics  equations, 
equations  (86)  through  (89)  .  The  plate  established  autorotation  within 
two  revolutions  as  depicted  in  Figure  21.  The  trace  of  the  numerically 
predicted  flight  path  began  at  the  point  where  the  card  was  "released" 
and  allowed  to  free  fall.  Experimental  data  of  a  plate  with  identical 
physical  characteristics  obtained  using  high  speed  photography  (Ref  58) 
is  shown  in  Figure  22.  The  experimental  flight  path  was  photographed 
after  autorotation  was  fully  established  and  thus  did  not  include  data 
for  the  launch  phase. 
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Trie  calculation  of  the  free  fal'-.nc,  autO‘'c  t;  tine,  flat  plate 
correctly  predicted  tne  complex  relationship  between  the  plcte's 
rotational  velocity  and  the  linear  velocity  of  its  center  ox  gravity. 
The  experimental  data  presented  in  Figure  22  illustrates  how  the  plate 


appears  to 

pivot  about 

its 

endpoints.  The 

magnitude  of 

the 

nondimensional 

1  inear 

velocity 

of  the  center  of 

gravity  and 

the 

nondimensional 

rate 

of 

rotation  were  changing 

periodical ly 

c  S 

illustrated  in 

Figures 

23 

and  24 

respectively.  The 

angle  of  attack 

i  r. 

these  figures  was  obtained  using  equation  (72).  The  relationship  between 
the  phasing  and  magnitude  of  the  angular  velocity  and  linear  velocity 
must  be  precise  to  present  the  appearance  that  the  plate  rotates 
alternately  about  its  endpoints.  This  complex  phenomenon  was  accurately 
predicted  by  the  numerical  calculations  indicating  the  procedure 
correctly  modeled  the  complex  flow  field  and  properly  related  the  flow 
field  to  the  three  degrees  of  freedom. 

The  plate  in  the  numerical  simulation  was  released  at  3.0  ft/sec 
and  developed  into  autorotation  with  a  mean  velocity  of  approximately 
3.01  ft/sec.  The  mean  velocity  of  the  experimental  plate  was 
approximately  3.75  ft/sec.  The  plate  in  the  numerical  simulation 
autorotated  at  a  mean  rotational  velocity  of  19.56  rad/sec  while  the 
experimental  plate  was  found  to  rotate  at  approximately  23.4  rad/sec. 

The  predicted  flight  path  of  the  plate  did  not  compare  well  with 
the  experimental  flight  path.  The  calculations  predicted  a  mean  flight 
path  angle  of  approximately  -0.20  radians  compared  to  -0.81  radians  in 
the  experimental  data.  This  discrepancy  in  flight  path  angle  may  have 
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peer.  causec  b;  znree  cirensional  effects.  The  experimental  plate 
have  produced  a  complex  three  dimensional  vortex  which  degraded  the  lift 
to  crac  ratio.  The  calculations  simulated  a  plate  with  an  infinite 
aspect  ratio  and  did  not  model  the  three  dimensional  effects.  The 
numerical  plate  should  therefore  have  a  higher  lift  to  drag  ratio  tnar 
the  experimental  plate  which  would  tend  towards  a  more  shallow  flight 
path  angle. 

The  calculated  aerodynamic  coefficients  produced  by  the  plate 
during  free  fall  autorotation  are  shown  in  Figure  25.  The  maximum  drag 
lags  the  maximum  lift  by  approximately  0.71  radians.  The  moment  is  in 
phase  with  the  lift.  In  contrast  to  the  plate  undergoing  forced 
rotation,  the  moment  is  not  in  phase  with  the  lift.  The  moment  being  in 
phase  with  the  lift  seems  to  be  a  necessary  condition  for  autorotation. 
The  moment  averaged  over  each  cycle  is  zero  which  is  required  in  fully 
established  autorotation.  The  experimentally  determined  coefficients 
for  a  plate  with  the  same  physical  characteristics  (Ref  58)  are  shown  in 
Figure  26. 

The  aerodynamic  coefficients  encountered  in  the  free  falling 
autorotaiion  calculations  were  considerably  larger  than  the  coefficients 
obtained  at  static  angles  of  attack.  A  comparison  between  the  static 
drag  polar  for  a  flat  plate  and  the  drag  polar  for  a  free  falling, 
autorctating,  flat  plate  are  shown  in  Figure  27.  The  large  differences 
between  the  static  and  dynamic  aerodynamic  coefficients  are  the  primary 
reason  that  past  attempts  to  calculate  autorotation  using  quasi-steady 
numerical  techniques  have  failed.  The  unsteady,  nonlinear  aerodynamic 
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c'-t  ar  essential  ingredient  ir.  cuto^ctatior.. 

Fur  the  first  time  since  initial  attempts  were  made  to  explain  the 
*>ee  falling,  autorotati ng ,  fiat  plate  in  1854,  the  complete  flow  fielc 
around  the  plate  is  known.  The  instantaneous  streamlines  and 
ccrresDcnding  lines  of  constant  vorticity  at  several  angles  of  attack 
arc  snown  in  Figure  28  and  29,  respectively.  Each  plot  has  been  rotated 
so  the  freestream  velocity  is  always  coming  from  the  right.  It  is 
difficult  to  see  the  vortices  being  shed  from  the  plate  in  Figure  28. 
However,  the  vorticity  contours  shown  in  Figure  29  illustrate  the 
asyrmetry  caused  by  the  rotation.  A.M.O.  Smith  realized  in  1953  that 
this  asymmetry  is  essential  for  autorotation. 

The  flight  paths  of  the  plate  released  with  different  initial 
conditions  were  also  found  to  compare  qualitatively  with  experiment. 
The  plate  was  released  similarly  to  the  previous  case,  however,  the 
initial  velocity  was  increased  to  6.0  ft/sec.  The  initial  flight  path 
was  slightly  different  as  shown  in  Figure  30  but  the  fully  oeveloped 

autorotation  was  identical  to  the  previous  case. 

In  another  launch  mode  the  plate  was  initially  moving  to  the  right 
with  a  velocity  of  6.0  ft/sec  at  a  glide  path  angle  of  -0.1  radians  and 
a  fixed  angle  of  attack  of  0.1  radians.  The  plate  was  "released"  and 
reversed  direction  as  shown  in  Figure  31.  The  plate  did  not  have 

sufficient  inertia  to  rotate  past  the  portion  of  the  flight  where  it 

experienced  a  retarding  moment.  The  plate  stopped  rotating,  reversed 
direction,  and  established  autorotation,  falling  in  the  opposite 
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di recti  or .  After  autorotetion  was  estabiis'-c,  t»  e  •‘"•g-t 
characteristics  of  the  plate  were  identical  to  tre  previous  free  fall 
autorotation  cases  only  the  plate  was  moving  in  the  opposite  direction. 
A  similar  flight  path  was  demonstrated  experimentally  b>  the  author, 
however,  it  was  not  documented  with  high  speed  photography.  Numerical 
experimentation  with  the  different  launch  modes  indicates  the  plate  will 
establish  the  same  autorotation  mode  independently  of  the  launch  mode. 
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SECTION  VII 


CONCLUSIONS  AND  RECOMMENDATIONS 


The  goal  of  this  research  was  to  develop  a  method  of  coupling  an 
unsteady,  nonlinear  aerodynamic  prediction  technique  with  the  flight 
mechanics  equations  and  to  show  the  method  was  a  valid  approach  to  use 
in  predicting  the  flight  path  of  a  free  falling,  autorotating ,  flat 
plate.  A  numerical  technique  which  combined  a  two  dimensional 
Navier-Stokes  aerodynamic  prediction  program  with  a  three  degrees  of 
freedom  trajectory  program  was  developed  and  used  to  predict  the  flight 
path  of  a  free  falling,  autorotating,  flat  plate.  The  overall  approach 
was  validated  by  comparing  the  numerically  predicted  flight  path  with 
the  experimental  flight  path. 

An  numerical  method  which  coupled  a  Navier-Stokes  equations  solver 
and  a  flight  mechanics  equations  solver  was  developed.  The  two 
dimensional  incompressible  Navier-Stokes  equations  were  solved  using  an 
implicit,  finite  difference,  numerical  solution  technique.  The  physical 
domain,  which  was  modeled  by  a  patched,  exponentially  stretched, 
Cartesian  grid,  was  transformed  intc  an  evenly  spaced  computational 
domain.  The  Navier-Stokes  equations  were  formulated  in  a  coordinate 
system  attached  to  the  plate  and  the  corresponding  velocities  were 
formulated  in  an  inertial  reference  system.  The  translation  and 
rotation  of  the  plate  were  modeled  by  imposing  them  as  a  boundary 
conditions  on  the  body.  The  resulting  system  of  equations  was  solved 
using  a  successive-over-relaxation  iteration  technique.  The  aerodynamic 
coefficients  were  calculated  by  summing  the  viscous  and  pressure  forces 
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around  the  plate. 


The  three  degrees  of  freedom  flight  mechanics  equations  were 
reduced  to  four  simultaneous  first  order  ordinary  differential 
equations.  The  aerodynamic  coefficients  were  entered  as  forcing 
functions  into  the  flight  mechanics  equations.  The  system  of  equations 
was  solved  using  an  open  Adams  method.  The  calculated  translation  and 
rotation  of  the  plate  were  used  as  new  boundary  conditions  in  the 
Navier-Stokes  equations  solver.  The  procedure  was  repeated  thousands  of 
times  to  generate  the  flight  path  of  the  plate. 

The  Navier-Stokes  equations  solver  was  used  to  predict  the  flow 
field  around  a  flat  plate  at  several  fixed  angles  of  attack.  The  flow 
field  was  unsteady  above  approximately  0.15  radians  angle  of  attack 
which  corresponded  to  a  lift  coefficient  of  approximately  0.44.  The 
maximum  static  lift  coefficient  for  a  flat  plate  was  found  to  be 
approximately  1.2. 

The  flow  fields  around  a  flat  plate  undergoing  forced  rotation  at 
three  different  rates  of  rotation  were  calculated.  The  maximum  lift 
coefficients  were  approximately  3.2,  5.8  and  8.5  for  nondimensional 
rotational  rates  of  1.0,  2.0,  and  4.0  respectively.  The  maximum  drag 
coefficient  and  the  maximum  moment  coefficient  lagged  the  maximum  lift 
coefficient. 

A  flat  plate  undergoing  pinned  autorotation  was  modeled  by  allowing 
one  degree  of  freedom  (rotation).  The  plate  was  released  at  a  small 
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ancle  cf  attack  and  established  autoroteti on  at  a  near  r.ondi mens i oral 
angular  velocity  cf  1.5S.  The  magnitude  of  the  aerodynamic  coefficients 
were  comparable  to  these  produced  by  the  plate  undergoing  forced 
rotation  at  an  angular  velocity  of  2.0.  The  maximum  drag  coefficient 
again  lagged  the  maximum  lift  coefficient,  however,  the  maximum  moment 
coefficient  was  in  phase  with  the  maximum  lift  coefficient. 

Finally  the  complete  set  of  equations  was  used  to  calculate  the 
flight  path  of  a  free  falling,  autorotating,  flat  plate.  The  flight 
path  predicted  by  the  numerical  method  compared  well  with  the 
experimental  flight  path.  The  maximum  lift  coefficient  produced  during 
fully  developed  autorotation  was  approximately  6.2  which  was  several 
times  larger  than  the  maximum  lift  coefficient  obtained  at  a  fixed  angle 
of  attack.  This  explains  why  previous  attempts  were  unsuccessful  in 
using  quasi-steady  aerodynamic  techniques  to  predict  autorotation.  The 
maximum  moment  was  in  phase  with  the  maximum  lift  which  seems  to 
indicate  the  moment  must  be  in  phase  with  the  lift'  to  support 
autorotation. 

Developing  and  demonstrating  a  procedure  to  solve  a  problem  which 
incorporates  the  essential  elements  of  a  spinning  aircraft  is  a  path 
finder.  This  path  finder  will  assist  in  the  future  analyses  of  the 
complete  aircraft  spin  problem.  In  addition,  this  research  effort  has 
advanced  the  state  of  the  art  in  several  other  areas: 

1)  Identifying  the  "mixed"  reference  system  for  future  application 
to  spin  computations  was  the  most  important  contribution  of  this  study. 
Previous  aerodynamic  calculations  for  rotating  bodies  were  performed  in 


coordinate  and  velocity  reference  systems  attached  fc  the  body.  Tr, i s 
approach  was  found  to  have  undesirable  characteristics  as  faster 
rotation  rates  were  introduced.  Hore  accurate  results  were  obtained 
when  the  coordinate  system  was  attached  to  the  rotating  plate  and  the 
corresponding  velocities  were  formulated  in  an  inertial  reference 
system. 

2)  Lugt  solved  the  two  dimensional  Ha vier-Stokes  equations  for  a 
thin  ellipsoid  rotating  at  a  constant,  predetermined  rate  of  rotation. 
However,  it  was  not  until  this  research  effort  that  the  rate  of  rotation 
w^s  allowed  to  vary  during  the  calculation.  The  rate  of  rotation  was 
calculated  from  the  aerodynamic  moment  which  resulted  in  true  pinned 
autorotation. 

3)  First  principle  aerodynamic  calculations  have  not  been  used  when 
solving  the  three  degrees  of  freedom  equations  of  motion.  A  method  was 
developed  to  couple  a  two  dimensional  Navier-Stokes  equations  solver 
with  a  three  degrees  of  freedom  equations  of  motion  solver.  To  the 
author's  knowledge,  this  was  the  first  time  an  unsteady,  nonlinear 
aerodynamic  prediction  technique  was  used  in  solving  the  three  degrees 
of  freedom  equations  of  motion.  The  overall  approach  was  verified  by 
comparing  the  numerical  results  with  experimental  data. 

4)  The  time  steps  used  in  finite  difference  solutions  of  the 
Navier-Stokes  equations  are  usually  constant  throughout  a  calculation  or 
are  changed  manually  by  the  programmer.  This  study  was  the  first  time 
the  time  step  in  a  Navier-Stokes  equations  solver  was  varied  during  a 
calculation  based  on  the  results  of  the  movement  of  a  body.  This 
significantly  reduced  the  computing  time  necessary  to  obtain  a  solution. 
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It  is  recommended  that  the  pursuit  of  i  numerical  solution  of  a 
complete  spinning  aircraft  be  continuec  ever,  though  the  final  ooal  is 
not  feasible  on  today's  computers.  Three  near  term  research  efforts  are 
recommended  to  continue  building  the  technology  base  necessary  to 
eventually  calculate  the  flight  path  of  a  spinning  aircraft: 

1)  The  cause  of  the  discrepancy  between  the  numerically  predicted 
and  experimental  flight  paths  should  be  determined.  The  numerical 
procedure  should  be  modified  to  improve  the  correlation  between  the  two 
fl ight  paths. 

2)  The  six  degrees  of  freedom  autorotation  of  a  three  dimensional 
object  should  be  calculated  using  linear  aerodynamics.  The  spinning 
maple  seed  offers  a  good  physical  example  of  such  a  problem.  The  flow 
field  around  the  maple  seed  is  probably  attached,  therefore  lending 
itself  to  analysis  by  linear  aerodynamics.  This  would  allow  the 
researcher  to  couple  the  two  systems  of  equations  without  expending  the 
large  amount  of  computer  resources  to  support  a  Navier-Stokes 
aerodynamic  solution. 

3)  Finally,  a  three  dimensional  Navier-Stokes  equations  solver 
should  be  coupled  to  a  six  degrees  of  freedom  equations  of  motion 
solver.  The  computational  method  should  be  used  to  calculate  the 
autorotation  of  a  relatively  simple  aerodynamic  shape  such  as  a  pure 
delta  wing.  This  problem  can  probably  be  solved  on  today's  computers. 
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Figure  1  Path  of  Free  Falling  Autorotafing  Flat  P'itt 


c)  Translation  and  Rotation 


Figure  7  Examples  of  Boundary  Conditions 
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Figure  10  Unsteady  Shedding  Vortex  Behind  Flat  Plate  at  a=0.4 
radians  and  Res5000 
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Figure  21  Calculated  Plate  Position  During  F 
Autorotation,  y  =-n/4,  V  '3.0,  0  =i 


Figure  23  Nondimonsi 


Figure  24  Nondimonslona 1  Rate  of  Rotation  of  a  Flat  Plate 
In  Free  Fall  Autorotation 
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Figure  27  Comparison  of  Static  and  Dynamic  Drag  Polars 


APPENDIX  k 


BODY  FIXED  ROTATING  SYSTEM 

The  Navier-Stokes  equations,  (14)  through  (16),  were  modified 
slightly  when  expressed  in  a  body  fixed  coordinate  system  undergoing 
arbitrary  pitching.  Arbitrary  translation  was  not  included  at  this 
point  (Refs  59  and  60).  The  body  fixed  flavier-Stokes  equations  are: 

ut  +  uux  +  vuy  =  -Px  +  h  <uxx  +  Uyy)  +  &  +  2©v  +  0y  (Al) 

vt  +  uvx  +  wy  =  -Py  +  ~  (Vxx  +  vyy)  +  0?  -  20U  -  0X  (A2) 

u  +  v  =  0  (A3) 

x  y 

The  far  field  boundary  conditions  for  velocity  and  pressure  were 
given  by  the  freestream  flow  field  in  the  inertial  reference  frame. 

These  conditions  are  equations  (A4)  through  (A6)  when  written  in  the 

body  fixed  coordinate  system. 

u  =  u  cos  (a  -  0)  +  0y 

a>  O  J 

V  =  Un  cos  (a  -  0)  -  0x 

P  =  P. 

The  no  change  boundary  condition,  which  allows  a  wake  to  pass  through 
the  boundary,  was  imposed  on  the  downstream  boundary  velocity  component. 
Numerical  experimentation  showed  that  imposing  freestream  boundary 
conditions  on  the  downstream  boundary,  instead  of  the  no  change 
condition,  had  little  effect  on  the  aerodynamic  characteristics  of  the 
plate.  The  downstream  outer  boundary  conditions  on  pressure  retained 
the  freestream  value.  In  the  fixed  system,  these  conditions  are  given 


(A4) 

(A5) 

(A6) 
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bv 


(A7) 

(A8) 

(a?; 


p  =  p. 


where  the  derivatives  are  normal  to  the  outer  boundary. 

The  flow  field  around  a  plate  undergoing  forced  rotation  was 
calculated  using  a  numerical  technique  similar  to  the  one  developed  by 
Hegna  (Ref  25).  The  solutions  were  acceptable  for  slow  rotational 
rates.  The  results  of  the  calculations  became  increasingly  worse  as  the 
rotational  rates  approached  a  reduced  frequency,  oc/U^k,  of 
approximately  one.  The  solution  remained  numerically  manageable  with 
the  evolution  of  time  even  though  the  pressures  near  the  outer  boundary 
approached  values  one  to  two  orders  of  magnitude  too  large. 

A  calculation  of  a  rapidly  rotating  flow  field  wfthout  a  body 
produced  similar  pressures  near  the  outer  boundary.  This  indicated  the 
problems  were  probably  caused  by  the  outer  region  and/or  the  boundary 
conditions.  The  same  calculation  was  performed  with  the  outer  boundary 
moved  farther  away  from  the  plate;  nearly  the  same  results  were 
obtained. 

Hodge  (Ref  41)  found  the  junction  between  the  upstream  Dirichlet 
and  downstream  Neumann  boundary  conditions  was  a  source  of  a  small 
error.  This  error  was  amplified  in  the  calculation  with  rotation.  The 
small  changes  in  the  velocity  components  necessary  to  satisfy  the 
boundary  conditions  and  continuity  at  the  intersection  of  the  upstream 
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end  downstream  boundary  conditions  became  unacceptably  large.  The  smell 
changes  produced  a  pressure  wave  which  propagatec  into  the  flow  fielc  at 
each  time  step.  In  the  steady  calculations,  the  flow  field  was  allowed 
to  stabilize.  Eventually  the  intersection  of  the  boundary  conditions 
was  a  source  of  only  a  small  error.  However,  in  the  rotating  cases  the 
source  of  error  was  continuous  and  increased  as  the  rate  of  rotation 
increased.  Averaging  and  smoothing,  which  considerably  complicated 
implementing  the  outer  boundary  conditions,  were  tried  as  methods  to 
alleviate  the  problem  but  to  no  avail. 

The  velocities  near  the  outer  boundaries  became  very  large  as  the 
rotation  rates  became  large  (Ref  61).  One  approach  to  alleviating 
numerical  problems  in  many  numerical  investigations  is  to  move  the  outer 
boundary  farther  away  from  the  body  (Ref  17).  In  this  formulation, 
moving  the  outer  boundary  increased  the  velocities  in  the  outer  region, 
thereby  compounding  the  problems.  The  grid  spacing  was  coarse  in  the 
radial  direction.  The  grid  spacing  was  equally  coarse  in  the  tangential 
direction  or  as  fine  as  the  minimum  grid  spacing  used  during  the 
problem.  The  exponential  grid  and  rotating  formulation  combined  large 
and  small  velocity  components  over  large  aspect  ratio  finite  difference 
grid  cells.  This  combination  occurs  at  the  location  of  the  largest 
velocities.  The  mismatch  in  grid  spacing  and  velocities  did  not  offer 
numerical  compatibility,  however,  it  did  not  directly  attribute  to 
numerical  difficulties. 

Hegna  (Ref  25)  was  successful  in  using  this  approach  to  calculate 
the  flow  field  around  slowly  rotating  airfoils  at  relatively  large 
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Reynolds  numbers.  Very  snail  time  steps  were  used  tc.  obtain  numerical 
stability  at  large  Reynolds  numbers  for  flow  fielcs  without  rotation. 
This  also  helped  stabilize  the  rotating  problem.  Usirc  a  time  step 
similar  in  magnitude  to  the  one  used  by  hegna  was  not  a  practical 
alternative  for  this  research  project.  The  computational  time  necessary 
to  do  several  rotations  of  the  plate  would  have  beer,  prohibitive. 

Finally,  a  Von  Neumann  stability  analysis  performed  at  certain 
locations  within  the  flow  field  showed  the  rotational  terms  could 
destabilize  equations  (Al)  and  (A2).  A  Fourier  representation  of  each 
variable  was  substituted  into  the  momentum  equations.  The  flow  field 
was  assumed  to  be  rotating  at  a  constant  velocity.  The  rotational  terms 
were  treated  as  source  terms  in  the  Navier-Stokes  equations  solver.  In 
the  stability  analysis,  however,  the  rotational  terms  were  formulated  as 
functions  of  u  and  v  by  assuming  that  u  and  v  were  of  the  same  order  of 
magnitude.  The  viscous  terms  were  assumed  to  be  small  in  the  flow  field 
away  from  the  plate.  The  stability  of  one  SOR  iteration  w as  found  to  be 
a  function  of  the  rate  of  rotation.  Above  some  rotational  velocity  the 
rotational  terms  will  destabilize  the  equations. 

The  rotating  reference  approach  was  abandoned  because  of  .  the 
difficulties  outlined  in  this  appendix.  A  "mixed"  system  was  adopted  as 
described  in  Section  III.B.  A  Von  Neumann  stability  analysis  showed  the 
"mixed"  system  had  better  stability  characteristics  than  the  rotating 
formulation.  A  plate  could  rotate  three  times  faster  in  the  "mixed" 
system  and  maintain  the  same  stability  characteristics  as  the  rotating 
formulation. 
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APPENDIX  E 


SOR  SWEEP  METHOD 


Calculating  the  flight  path  cf  a  free  falling  flat  plate  was  truely 
time  dependent  which  required  more  attention  than  normal  to  be  paid  to 
the  successive  over  relaxation  (SOR)  sweep  direction.  The  sweeping 
direction  was  important  because  it  effectively  produced  a  numerical 
propagation  speed.  This  speed  was  not  constant  throughout  the  flow 
field  because  of  uneven  physical  step  sizes.  Alternating  the  sweeping 
direction  allowed  the  numerical  propogation  speed  to  be  similar  in  all 
directions.  An  even  number  of  iterations  were  used  at  each  time  step  to 
avoid  favoring  one  direction.  An  “infinite"  numerical  propagation  speed 
more  closely  approximated  the  infinite  speed  of  sound  of  incompressible 
flow. 


The  concept  of  numerical  propagation  speed  is  illustrated  in 
Figure  Bl.  A  disturbance  was  imposed  at  point  5  in  the  flow  field  shown 
in  Figure  Bla.  One  iteration  of  a  calculation  used  a  left  to  right 
sweep  which  only  allowed  the  influence  of  the  disturbance  to  be  felt  at 
point  4  as  illustrated  in  Figure  Bib.  A  second  iteration  allowed  the 
disturbance  to  be  felt  at  point  3  and  so  forth.  Four  iterations  were 
required  to  propagate  the  disturbance  from  point  5  over  to  point  1  using 
a  left  to  right  sweep.  Figure  Blc  illustrates  that  a  calculation 
performed  with  a  right  to  left  sweep  allowed  the  disturbance  to 
propogate  throughout  the  complete  flow  field  in  one  iteration. 

The  most  ideal  SOR  sweep  mode  would  have  been  alternating  the  sweep 
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a)  Initial  Conditions 

1  1  i 


b)  Sweeping  to  the  Right 


c)  Sweeping  to  the  Left 


Figure  B1  Concept  of  Numerical  Propagation  Speed 
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outward  frorr.  the  plate  ir  ell  directions  and  inward  from  the  cuter 
boundary  in  ell  directions.  This  mooe  however,  would  have  been 
difficult  and  inefficient  to  implement  in  a  computer  program.  A  more 
efficient  sweep  mode  was  used.  During  the  first  iteration,  the  sweep 
started  at  the  middle  of  the  computational  domain,  then  swept  up  and 
down  from  the  plate  to  the  outer  boundaries.  The  starting  points  were 
moved  to  the  left  and  right  one  grid  space  and  swept  up  and  down  from 
the  plate  to  the  outer  boundaries.  This  process  was  repeated  until  the 
complete  domain  was  covered.  The  second  iteration  began  at  the  corners 
of  the  outer  boundary.  The  calculations  were  then  swept  to  the  center 
of  the  computational  grid  then  moved  back  to  the  outer  boundaries.  The 
starting  points  were  moved  in  one  grid  space  and  the  process  repeated 
until  the  complete  domain  was  covered. 
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REFERENCE  SYSTEMS  FCf.  CAT A  PRESENTATION 


Selecting  a  reference  frame  in  which  to  present  unsteady  motion  is 
difficult.  Steady  motions  can  usually  be  presented  in  a  body  fixed 
reference  system  in  which  the  flow  becomes  steady.  In  this  steady  state 
system,  streamlines,  pathlines  and  streaklines  coincide.  Unfortunately, 
in  unsteady  flows,  no  preferred  reference  frame  exists,  therefore,  the 
reference  frame  was  selected  on  the  basis  of  other  criteria. 
Streamlines,  pathlines  and  streaklines  generally  do  not  coincide  in 
unsteady  flow  fields  (Ref  62). 

The  simplest  case  of  a  rotating  body  in  two  dimensions  is  that  of  a 
flat  plate  rotating  with  constant  angular  velocity  in  a  fluid  at  rest  at 
infinity.  Figure  Cla  shows  the  approximate  instantaneous  streamlines  of 
a  potential  flow  around  a  rotating  plate  (Ref  63).  This  flow  is 
unsteady,  but  by  superimposing  the  rotation  of  the  plate,  the  flow 
became  steady  as  shown  in  Figure  Clb. 

The  flow  becomes  more  complicated  when  the  plate  rotates  with 

constant  rotation,  ft  ,  in  a  uniform  freestream,  Va  ,  because  no  steady 

state  exists.  Four  different  frames  are  distinguishable  (Ref  62).  When 
*  * 

U  and  n  are  the  translational  velocity  and  rotation  of  the  plate 

relative  to  the  reference  frame,  these  four  different  frames  are: 

*  ★ 

1)  U  *0,  0  =0;  the  frame  is  fixed  to  the  plate. 

*  * 

2)  U  ft  =0;  the  frame  does  not  rotate  relative  to  the  plate; 

but  translates  relative  to  the  plate. 
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★  *■ 

3)  L'  =C,  L  the  frame  is  fix.ec  tc  the  plate  v.-ith  !-egerd  to 
translation;  but  rotates  relative  to  the  plate. 

*  it 

4)  U  =V  ,  !7  =fi;  the  frames  is  in  translational  and  rotational 

7  co 

motion  relative  to  the  plate. 

The  streamlines  around  a  free  falling,  autorotating,  flat  plate  at 
a=l£.34  radians  angle  of  attack  for  the  four  different  cases  ere 
illustrated  in  Figure  C2.  Since  the  streamlines  are  not  invariant,  all 
patterns  differ  from  each  other.  Frame  1  is  advantageous  since  the  body 
contour  was  a  streamline  and  the  flow  field  near  the  plate  could  be 
examined.  Frame  3  is  appealing  because  of  its  analogy  to  instantaneous 
streamline  photos  of  a  rotating  plate  in  a  wind  tunnel.  Frames  2  and  4 
are  of  less  interest  in  this  study. 
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d)  Frame  4 


Figure  C2  (cont.) 


APPENDIX  i 


ADAK-BASHFORTH  FORMULAE 


The  Adam's  multistep  formulae  have  been  used  in  efficient  computer 
routines  to  obtain  the  solution  of  ordinary  differential  equations.  The 
Adam's  formulae  fall  into  two  general  categories:  open  formulae  and 
closed  formulae.  The  Adam's  open  formulae,  also  known  as  the 
Adams -Bashforth  formulae,  were  the  only  formulae  used  to  solve  ordinary 
differential  equations  in  this  study. 

The  open  formulae  were  derived  by  performing  a  forward  Taylor 
series  expansion  about  an  arbitrary  value  of  time.  The  derivatives  in 
the  expansion  were  replaced  by  various  numbers  of  backward  differences 
depending  on  the  order  of  the  formula.  The  coefficients  of  different 
order  Adam-Bashforth  formulae  are  shown  in  Table  D.I. 

These  formulae  are  called  open  formulae  since  the  next  value  in 
time  of  the  function  could  be  solved  for  explicitly  in  terms  of  the 
present  and  past  values  of  the  function.  However,  during  start-up  of 
the  second  order  and  higher  order  methods,  the  previous  values  of  the 
function  were  not  known.  The  previous  values  of  the  function  were 
obtained  by  linearly  extrapolating  back  in  time.  These  extrapolated 
values  made  the  Adam's  formulae  self  starting  regardless  of  the  order. 
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Table  D.I.  Coefficients  of  the  Adams  Formulae 
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